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Summary 

An investigation is conducted of several numerical schemes for use in the computation of two- 
dimensional, spatially evolving, laminar, variable- do ns it y compressible shear layers. Schemes with various 
temporal accuracies and arbitrary spatial accuracy for both inviscid and viscous terms are presented 
and analyzed. All integration schemes use explicit or compact finite- difference derivative operators. 
Three classes of schemes are considered: an extension of MacOormack s original second-order temporally 
accurate method, a third-order temporally accurate variant of the coupled space-time schemes proposed 
by Rusanov and by Kutler, Lomax, and Wanning (RKIAV), and third- and fourth-order Runge-Kutta 
(RK) schemes. The RKLW scheme offers the simplicity and robustness of the MacCormack schemes and 
gives the stability domain and the nonlinear third-order temporal accuracy of the Runge-Kutta method. 
In each scheme, stability and formal accuracy is considered for the interior operators on the convection- 
diffusion equation 1 7 + al r j- = o r l . for which a and o r are constant. Both spatial and temporal 
accuracies are verified by the equation Ut — [ft(.r)r>].r as well as I f + bj — 0. Numerical houndary 
treatments of various orders of accuracy are chosen and evaluated for asymptotic stability. Formally 
accurate boundary conditions are derived for explicit sixth-order: pent adiagonal sixth-order; and explicit, 
tridiagonal, and pent adiagonal eighth-order central-difference operators when used in conjunction with 
Runge-Kutta integrators. Damping of high wave-number, nonphysical data is accomplished for all schemes 
with explicit filters, derived to tenth order on the boundaries and twentiet h order in the interior. Several 
schemes are used to compute variable-density compressible shear layers, where regions of large gradients 
of flow-field variables arise near and away from the shear-layer centerline. Results indicate that in the 
present simulations, the effects of differences in temporal and spatial accuracies between the schemes are 
less important than the filtering effects. Extended MacCormack schemes are robust but inefficient because 
of restrictive Courant -If iedrichs-Levy (OFL) limits. I he third-order temporally accurate RKIAV schemes 
are less dissipative but have shorter run times. The Runge-Kutta integrators did not have sufficient 
dissipation t o be useful candidates for the computation of variable-density compressible shear layers at the 
levels of resolution used in t he current work. 

Introduction 

The numerical simulation of spatially evolving, compressible shear layers has become popular as a tool 
to understand the mixing mechanisms involved in supersonic combustion. Simulations may involve not only 
the effect of compressibility but also the presence of large gradients in density caused by disparate-mass gas 
mixtures or large temperature gradients that arise from exothermic chemical reactions. In disparate- mass 
gas mixtures, the Schmidt and Lewis numbers are non unity- usually greater than 1 in one stream and 
less than 1 in the other. Self-similar solutions to the laminar shear layer suggest that this nonunitv gives 
rise to different profiles for species and temperatures relat ive to the velocity profile. In hydrogen-nitrogen 
mixing layers, vorticity occurs predominately in the low-density stream. (See ref. 1.) Ibis phenomenon is 
experimentally observed in turbulent, disparate- mass supersonic shear layers. ( See ref. 2.) 

Computation of compressible shear layers has been largely confined to gas streams that are uniform in 
composition; also, the specific numerical method chosen has varied considerably. Soetrisno et al. (ref. .5) 
use a second-order accurate, finite- difference, total variation diminishing (TVD) scheme coupled with a 
second-order accurate Runge-Kutta method to study two-dimensional, temporally evolving, inviscid shear 
layers. Yamamoto and Daiguji (ref. 4) use either a fifth-order upwind TVD or a fourth-order monotonic 
upwind -centered scheme for conservation laws (Ml'SCL) TVD scheme, coupled with a ( ran k-Ni colson 
time integrator. Shu et al. (ref. 5) use various order, essentially nonoscillatory (ENO) finite-difference 
schemes, as well as compact central -difference stencils and a third-order low-storage' Runge-Kutta. method 
on a t hrtv-dimensioiial shear layer. The ENO schemes are particularly useful for flows in which steep 
gradients are present. Grin stein and Kailasanath ( ref. 6) use a flux-corrected transport (F(T) algorithm 
to investigate three-dimensional and chemical-reaction effects. Another method that has been us*' fill in the 
simulation of compressible flows is upwind biased differencing. Raiand Moin (ref. i ) use mildly dissipative 
fifth-order upwind differences on inviscid terms and fourth-order differences on viscous terms together 



with an implicit time integration to simulate transition and turbulence in supersonic boundary layers. 
In a different approach to adding dissipation, Mukunda ct al. (ref. 8) use the (2-4) scheme proposed by 
( lOttlieb and Turkel (ref. 9), which is second-order accurate in time and fourth-order accurate in space, as 
well as the compact (2-4) version of the Mac Corn lack (ref. 10) method developed by Carpenter (ref. 11) 
to study spatially evolving compressible shear layers. Lele (ref. 12) chooses a sixth-order compact, central- 
difference stencil for viscous and in vise id terms and a third-order low-storage Runge-Kutta method to 
calculate temporally and spatially evolving, two-dimensional compressible shear layers. Dissipation is 
added by the use of implicit filters. (See ref. 13.) In a combination of compact finite- difference and Fourier 
spectra] methods, Sand ham and Reynolds (ref. 14) investigate the transition of a compressible shear layer; 
Guillard, Male, and Peyret (ref. 15) use a fully spectral scheme. 

Computations may be divided into two broad categories - spatial and temporal simulations. Temporal 
simulations allow the use of periodic numerical and physical boundary conditions in the streamwise 
direction, which greatly simplifies the computations. Unfortunately, they are an idealization of real shear 
layers. Spatial simulations require specification of both physical and numerical boundary conditions. 
Recently, Carpenter, (lottlieb, and Abarbanel (ref. 16) have determined numerical boundary treatments 
that preserve the accuracy of compact, tridiagonal sixth-order interior schemes on the model hyperbolic 
equation l : t + aUr — 0. These treatments are also asymptotically stable with respect to time. Prior to 
the work of Carpenter, researchers using the sixth-order tridiagonal stencil for the interior scheme closed 
it at the boundaries so that the formal accuracy of the overall method was reduced. 

Another relevant issue in the computation of compressible shear layers is numerical dissipation. In 
simulations where not all the relevant length and/or time scales of the problem are being resolved, 
dissipation must be added to ensure computational stability. Some numerical dissipation is desirable 
to remove spurious high-frequency information regardless of whether second-order derivatives are taken 
once with a second- order derivative operator or twice with a first-order derivative operator. The source of 
this high-frequency information may be intrinsic instability in the scheme, the misspecificat ion of physical 
houndary conditions, the “odd-even" decoupling between grid points, or insufficient resolution (temporal 
and spatial). To address this problem, some researchers have resorted to implicit (ref. 13) and explicit 
filters (refs. 17, 18, and 19). In the present work we use explicit filters. 

The goal of the present study is to generate families of schemes with arbitrarily high spatial accuracy 
for both viscous and inviscid terms, coupled with explicit time integrations from second to fourth order. 
Schemes are applied to highly resolved, spatially evolving, compressible shear-layer calculations devoid 
of discontinuities. The accuracy, stability, and robustness of the schemes are considered with particular 
attention to compressible, variable- density, noureacting flows. Analyses of both compact, and explicit 
interior schemes have been provided , as well as a variety of choices for boundary closures and explicit, filters. 
Stabilit y is considered not only through Von Neumann analysis but also through mat rix a nalysis of various 
boundary and interior treatments. The schemes examined are the following: extended MacCor mack -type 
schemes, anew variant of the schemes presented by Rusanov (ref. 20) and by Kutler, Lomax, and Warming 
(refs. 21, 22. and 23) (RKLW), and Runge-Kutta (RK) schemes. 


Numerical Method 


The governing equations are solved in conservative form with the SPAR.K2D (ref. 2-1 ) code and may be 
written as 


tfU 0F(U) 
dt dx + 


d>G(U) 

dy 


( 1 ) 


2 




fm 

pUU — <T rr 
pur — <Ty iV 

~ - (TfyV 

puYj + puiYi 


pr “I U 

pur - 0 

G = pVV — (T yy H= 0 

(pt 0 — <Tf fy)r — <T//.r W “h ([// 0 

pv V / T- pvj Y; _ _ 0 J > 

p is the density: ?/ is the streainwi.se velocity; r is the transverse velocity : (r.j n is the Newtonian stress 
tensor: f o is the total internal energy: q n is the heat, flux vector; V; is the species mass fraction; and it; 
and vj are st. teamwise and transverse components of the diffusion velocity, respectively. Homan indicies 
(e.g.. y) correspond to the species index, whereas Greek ones (e.g., o ) correspond to spatial indicies. 
Throughout this text, the inviscid derivative operators are those used to differentiate F and G. whereas 
viscous derivat ives are t hose used to generate de rival ive terms in the expressions for the stress tensor, heat 
flux vector, and diffusion velocity. 

In the finite- difference schemes considered here, for constant grid spacing A.r. tin' spatial derivative' of 
a function / (/' = /, r ). is given in matrix form as 

p f, = 4" q/ 

Ax 
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To avoid the increased operation count necessary to invert large bandwidth matrices, the bandwidth of the' 
matrix P is not considered operationally to be large' r than pent adiagonal . De term in at ion of t he specific 
cent ered-dilfere nee stencil is accomplished by writ ing 
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where the coefficients a, • • • form lilt' matrix P; o, b, ■ ■ ■ form the mat rix Q; and fj and /j are the values 
of some function and its derivative at grid point i, respectively. By defining the Fourier transform and the 
inverse Fourier transform of the discrete function value /,„ with (ref. !2;>) 
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where £ is the Fourier dual variable; if <F(£) is the approximation of the derivative of / in Fourier space or 
the Fourier image of / , /,„ is given ill Fourier space as 
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The finite-difference stencil given in equation (4) becomes 
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Consequently, 


/[2a sin(£) + 26 sin(2£) + 2c sin(3£) + 
[1 + 2o cos(<£) 4- 2/J cos(2£) + ■ • •] 


With the spectral representation of f(x) written as 
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it can he seen that / («r) in Fourier space has the form 
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If the Fourier image of the finite-difference derivative operator is expanded in a Taylor series in C coefficients 
of the stencil can be chosen to approximate the spectral derivative to some desired accuracy (i.e.. ^ % /£). 

For ail arbitrarily skewed stencil, the stencil and its Fourier image are given bv 
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where the subscripts L and R are used to denote left and right. 
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Filially, for predictor-corrector dissipative schemes {ref. 1(5). we have 
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The | > re diet or-comv to r stencils become centered stencils for B — D — f — G = 0. All second derivatives 
are taken effectively by successively applying a first derivative operator twice. A consequence of this 
application is that the wave number £ = n becomes neutrally stable for central difference schemes and 
may cause a. loss of stability on nonlinear problems. This wave number is sometimes referred to as the ~jr' 
mode. 

When predictor- correct or finite- difference schemes are used, viscous derivatives are calculated with 
explicit stencils (i.e., P = I. where I is the identity matrix, and Q = A). An explicit stencil of 
(N — l)th-order accuracy is used for the evaluation of viscous terms in the schemes when the derivatives 
of F and G are calculated to A- til-order accuracy. Runge-K ut ta schemes use the same derivative operator 
for both viscous and inviscid derivatives. Further discussion of the derivation of the stencils is contained 
in appendix A. 


Extended MacCormack Schemes 

in 19(59, MacCormack (ref. 10) introduced a two-stage numerical scheme for compressible flows with 
a predictor stage followed by a corrector stage, the scheme is second-order accurate in both space and 
time and is widely applicable, in part, because' of its simplicity and robustness. Details of the method 
can be found in many places (refs. 2(5 to 29). Attempts made by Gottlieb and lurk el (ref. 9) to improve' 
the method increased the inviscid spatial accuracy to fourth order. This scheme has been popular among 
researchers involved with highly resolved flow fields. (See refs. .10, .51, .52, and 8.) ( arpenter (ref. 11) 
further modified this scheme by using a compact fourth-order inviscid stencil with a third-order upwind 
viscous stencil. 1 he scheme was slight ly more accurate than the Gott liele 1 u rkel scheme. Bayliss (ref. .M) 
extended the GottlieleTurkel scheme to sixth-order accuracy for the inviscid terms. 

Extended MacCormack schemes take the original MacCormack scheme to arbitrary spatial accuracy in 
both the inviscid and viscous terms. These schemes are obtained by using the skewed stencils (eq. (12)) 
to generate viscous terms in the vectors F and G. and the predictor-corrector stencils (eq. (M)) are used 
to evaluate the derivatives of F and G. Symbolically, the schemes may be represented for the equation 


as 
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wliere A + and A arc the forward and backward difference operators and A = is the C'ourant- 
Friedriclis-Levy (CFL) number. 


The stability of the extended MacOormack scliein*'s may be conveniently analyzed in Fourier space with 
conventional Von Neumann analysis on the convection-diffusion equation V, + aU r = a r U rx with a and 
«r constants. If <F and — 'F* (where ty* is the complex conjugate of ♦) are defined as the Fourier images 
of A + and A - , <F v and — are defined as the Fourier images of A^ - and A r , the viscous derivatives; V 
is defined as the Fourier transform of U": \' = q$L-, A,, = ApH is the viscous CFL number or diffusion 
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Ls the amplification factor, then these schemes can be written as 
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The amplification factor G represents the magnitude of the amplification of a given frequency when 
the solution is advanced one time step. Use of the letter G here should not be confused with its use in 
defining the dissipative stencils. To determine the maximum CFL number, G may be analyzed for the 
interior scheme or the amplification matrix G may be considered for the full scheme with boundary points. 
For — 7T < £ < 7T, the magnitude of G must never exceed 1, or the spectral radius of G must be less than 
or equal to 1; this is required for stability. For the convection-diffusion equation U( + aUx — j j, the 
amplification factor and matrix are written as 


\G\ = 


“[!+(!+ A'** - A r **tf r ) (1 - A'* - A r #**)] 


(19) 


G = i[I+ (I-VA- + A.A-A+) (I - A'A+ + A r A+A")] (20) 

where A+ and A - are the forward and backward inviscid matrix derivative operators and Ajt* and A~ 
are the forward and backward matrix operators for the viscous derivatives. 

For consistency, the explicit coefficients must sum to zero as follows: 


(F - E) + (D - C) + (B - A ) + <7+ ( B + A) + (D + C) + (F + E) = 0 (21) 

or G = — 2(B + D + F). Values of B, D , and F may be selected on the basis of their effect on the 
dispersion, dissipation, and CFL numbers of the scheme. We have chosen to maximize the inviscid CFL 
number and retain the largest viscous CFL limit possible. The coefficients A, C\ £\ n, and /? are chosen 
to satisfy the accuracy requirement of the stencil. 

Table 1 lists the coefficients of the extended Mac. Cor mack predict or -cor rector stencils and the maximum 
C FL number Amax for the inviscid problem (ar = 0) in the absence of boundaries. The letters E. T, and 
P indicate that the matrix P is either diagonal/explicit, tridiagonal, or pentadiagonal, respectively. The 
notation (2-6E) should be interpreted as second-order temporal accuracy with sixth-order spatial accuracy 
for both the inviscid and viscous terms; the letter E indicates that the inviscid derivative operator is 
explicit. Some confusion may arise because many schemes found in the literature do not treat viscous 
terms at all and others do not ret ain the stated inviscid accuracy on viscous terms. 
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Figure 1 presents the stability boundary of these schemes as a function of the viscous and inviscid (T L 
numbers, again in the absence of boundaries. Regions in I lie lower left portion of the figure represent the 
stable region, whereas regions in the upper right portion are unstable. While optimizing these schemes 
it was noticed that contours of the stability boundary can be dramatically altered by different choices of 
/y, I ), and F. Optimizing was done by simply scanning parameter space for combinations with desirable 
stability characteristics. Each scheme can be considered as opt imized, although a o-percent improvement 
may be possible. Care should be exercised in optimizing the “flipping" parameters B, D, and F because 
many of the combinations share the maximum inviscid Cl L limit of the scheme yet very few of this subset 
have a boundary G = 1 that does not intersect the origin. The (2-8T) scheme was found to have no values 
of D . and F for which the boundary G - 1 did not intersect the origin. In many flow fields of interest . 
the local viscous and inviscid CFL numbers are likely to lie outside the stability domain in the (2-8T ) 
scheme. Explicit schemes possess significantly larger stability domains than their compact counterparts 
because of the increased truncation error of the explicit derivative operator. Thus, the CFL limits for 
compact schemes are more severe than those for multidimensional schemes. (See refs. 10 and lb) 

Note that the (2-4E) scheme differs from that proposed by Cottlieb and Turkel (ref. 9) where B — — A 
and I) = —C. If used with t he explicit third-order accurate viscous derivat ive, the Cott lieb-T urkel scheme 
has a viscous CFL limit of zero as CFL — 0. MacCormack's original scheme. (2-2E), is included in table 1 
and figure 1 for completeness. The explicit skewed viscous stencils ( /^ = o j = or = T ^ = 0) are given 
in t able 2 . 

RKLW Schemes 

Rusanov (ref. 20) derived a finite-difference scheme for nonlinear hyperbolic systems that, was uniformly 
third-order accurate in space and time. This scheme was considered for use in the computation of 
discontinuous solutions. T hree spatial difference operators were used in its construction mean value, 
difference, and identity. These operators were combined with a three-stage, third-order Runge-kutta 
method. Later, Burstein and Mir in (ref. 34) derived a similar method. Because function evaluations 
needed to be made on a staggered mesh, Kutler, Lomax, and Warming (refs. 21, 22, and 2 3) adapted 
Rusanov's scheme by replacing the first two stages with MacCbrmack s scheme. Hereinafter, this scheme 
is referred to as RKLW . This adaptation made the programming logic simpler and facilitated the inclusion 
of a source term and the extension to multidimensions. Various investigators have applied this scheme 
to both high-speed flow (ref. 35) and to meteorological flows (refs. TO. 37, and 38). Further discussion of 
the RKLW schemes can be found in the textbook by Anderson. TannehilL and Piet die r (ref. 26) and in 
two papers by Yanenko et al. (refs. 28 and 29). Attempts to proceed to uniformly fourth-order schemes 
for hyperbolic equations (refs. 39. 40. and 41) have been successful, but have not been used extensively, 
probably because' of their enormous complexity. 

The proposed RKLW method is a generalization of the third-order predictor-corrector format of Kutler, 
Lomax, and Warming (refs. 21 ,22, and 23) to arbitrary spatial accuracy in both viscous and inviscid terms 
within the temporally third-order Runge-Kutta (RIv) accuracy constraints. The implementation of this 
scheme for the equation I f + I\ v — 0 is solved numerically as 


r; = r/' - ;h \ a+f" 
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rf = (i-^)r/ , + /vr 
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where A r is the central -difference operator and the forward and backward differencing operators are the 
same as those' used in the extended MacCormack schemes. The values of /? j and ii-> are 1 and 1/4, 
respectively. The two degrees of freedom inherent in the general three-stage, third-order Runge-Kutta 
formulation (ref. 42) are used to accommodate symmetric (o;ji = o 32 in eq. (2,1)), predictor-corrector 
spatial differencing. 

The traditional Runge-Kutta scheme may be represented for the equation U t = —F* — = -fU* 

as 

U* - U" -o 2 l Xf n A l V n 
V ** = U n - a%\\f" A.\U' 1 - « 32 A f*A 2 F* 

,n>+\ _ v „ _ bl xf" Al U n - b- 2 \f*A 2 U* - IriXf** A-iU** 



where the subscript associated with the matrix operator A represents the finite difference operator used 
on the relevant stage —forward , backward, or centered. 

Runge-Kutta schemes are often described in terms of the Butcher array. (See ref. 42.) The Butcher 
array for the present scheme is given as 
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Hence, — <*21 anc ^ $2 ~ a 3l — <y 3*2- Symbols (letters) used for terms contained in the Butcher 
array should not be confused with those symbols (letters) involved in the definition of the finite- difference 
coefficients. Values of c; correspond to the time at which the ith stage is evaluated, that is, zero being the 
nth stage and one being the ( n -f 1 )th stage. 


The stability of the scheme is considered in Fourier space in the linear model equation Ut-\- aU* — ftvUr j 
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where is the Fourier image of the central-difference stencil. For the final stage of this method, the 
dissipative difference operator is used as a central- difference operator by setting B = D — F = 0. The 
optimum values of B, D, and F are different from those of the extended MacCormack schemes. They 
have been chosen to maximize the size of the stability envelope under the constraint that |(7| < 1 or that 
the spectral radius of G is less than or equal to 1. where 
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(2(5) 


a = i + |(i - v* r + A r * r # r ) 

x { 1 + ;h [( 1 + - /*i A,. ¥**,.) ( 1 - ,^iA> - ,i|A r 'M»J) - l] } 

G = il+|(I- A'A r + A l ,A , 'A r ) 

x { I + ,i 2 [(I - ;^A'A _ + rfiA,.A“ A+) (I - ;*,A'A + + A t A + A“) - i] } (27) 

Table *5 summarizes the RKLW schemes. Viscous stencils are the same as those used in t lie extended 
Mac( or mack schemes. (See table 2.) Figure 2 presents the stability boundary of these schemes in the 
absence of boundaries. As for the extended MacCormack schemes, explicit schemes have less restrict ive 
CFL bounds than the compact schemes with the same spatial order of accuracy, and increasing the order 
of accuracy of the spatial operator reduces the stability domain. The (3-2 F) scheme is the analog of the 
original MacCormack scheme with the Rnnge-Kutta integrator. By simply setting B — D — F = G — 0, 
the RKLW scheme becomes the third-order Runge-K lit t a scheme. Both RKLW and its corresponding 
third-order R unge-Kutt a. method require three storage locations (3M ) as opposed to the low-storage (2M) 
method proposed by Williamson (ref. 13). 

Runge-Kutta Schemes 

A significant portion of direct numerical simulations (DNS) and well-resolved model-free simulations 
of compressible flows have used a Rnnge-Kutta method. Unlike the extended MacCormack and RKLW 
schemes, space and time are not coupled in the numerical method. A commonly used Runge-Kutta method 
is the three-stage, third-order, low-storage scheme (ref. 43). Combined with a sixth-order, compact central- 
difference operator, this method has been used m the simulation of compressible shear layers (ref. 12), 
supersonic boundary -layer transition (ref. 44), and compressible isotropic turbulence (ref. 45). 

In a more traditional approach, common variants of the third- and fourth-order Runge-Kutta time- 
integration schemes are combined with explicit and compact differencing. The Butcher array for these is 
given by 



for third-order temporal accuracy and by 
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for Point li-ordor temporal accuracy. Tlie application of these schemes to the equation Ff + /’, = () is 


Ff = Uj - a-2\\y‘F" \ 

Ff* = Ff - ay 2 \A<F* l (30) 

f '" +1 = Ff + {biA'F" + b-,A r Ff + b$A r F**) J 

and 

Ff = Ff - a. n \A' F" 

Ff* = l’" - a :ii \A‘F* 

> (31) 

Ff = Ff - a^XA'F** 

Ff +l = Ff + (hA'Ff + b 2 A c Ff + b$A c F** + b 4 A r Ff**) , 

Analysis of the stability of the Runge-Kutta schemes can be done again with the equation 
Ff + aFj. = a r Fj.j.. The amplification polynomial (ref. 42) for the linear problem is given by 


G = M X> Z + E h ' c ' ) Z ' 2 ~ E Z ‘+ I E bowk ! Z 1 + ••• 


(••Vi) 


w'= 1 


FJ = 1 




where Z may be eit. her (A'4» r - ) or (X f A ( — X V A C A ( ), depending on whether one is interested in 

amplification factor or matrix, and ms is the number of stages in the Runge-Kutta scheme. The schemes 
are summarized in table 4. Figures 3 and 4 show the stability boundary of the third- and fourth-order 
Runge-Kutta /centered- difference schemes, determined by the amplification factor. Stability appears to 
be significantly augmented by going to fourth-order accuracy. Figure 3 shows a characteristic of tliird- 
and fifth-order Runge-Kutta. formulas: a tendency for the stability domain to become small as the viscous 
(TL — * 0. These stability domains are independent of which of the two free- parameter families of three- 
stage, third-order and four-stage, fourth-order Runge-Kutta schemes are chosen. All Runge-Kutta schemes 
considered use centered stencils to evaluate viscous derivatives. A brief discussion of low-storage Runge- 
Kutta schemes is contained in appendix B. The stability of Runge-Kutta schemes applied to the Navier- 
Stokes equations has been considered by Sow a (ref, 46) for second-order centered spatial derivatives and 
Runge-Kutta coefficients in which all (i}j = 0, except when i — j - 1-1. Temporal accuracy of the Runge- 
Kutta schemes was verified in the representative (4-8P) case. The linear equation Vf + l 7 x — 0 was 
solved for various (_TL numbers with a sinusoidal initial and boundary condition and 75 grid points per 
wavelength. As table 5 shows, fourth-order temporal accuracy was recovered. Table 6 contains the error 
when Vf + Uj- — 0 was solved at a OFL number of 0.01 on various grids to ensure that, spatial error 
dominated the total error. As can be seen, eighth-order spatial accuracy was recovered. In each of these 
cases, machine precision becomes a factor at liigli resolution. 


Formal Accuracy 

To determine the formal spatial and temporal accuracy of the interior schemes used in this study, the 
amplification factor and the linear equation Vf + aU x — <\ v V xx are used to derive the modified equation. 
(See refs. 26, 27, 47, and 48.) 


The exact solution to the convection- diffusion equation can be solved with the continuous Fourier 
transform as 


or 


j ^ I i 4“ iout ( T o l ^ t 


I f -f- iauV T o pu/ l — 0 


0 


(33) 

(34) 
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Tins equation may be integrated as 


In / 


where t u is the nth time step, to give 


/” + r 
t" 


— — ^/V/wV + o r a/*^ / |J/j 


+r 


r 


l n +T 


= exp + o r ^^ r f : |j„ 


At r — (A/), the value of the amplification factor is then 

f'Vxact — ex P + o r u;^ (A/) 

or. in terms of £, 




'exact — ex P — ^ A £ + A r £^ 


where A = ^77 and Ar = 7 -^,-, . 

- Ai (A.r)- 

r fhe error of the numerical scheme in Fourier space can now he written as (refs. 49 and 25) 


A/ 


1 {1 ( ^scheme ) 1 n ( G^xarl ) 


= v ,hm‘ 

k= 0 


A/ 


A- 0* 


Replacing all occurrences of £ h with its transform ( — /A.r )* ~~[.d m ) yields 


111 (f»«clieiiK‘) (G exact) 


y sdie iiK 

A / 


d k 


A/ 




n-=o 


A symbolic manipulator may now he used to expand In (i in a ]>ower series and to solve for A j.. 

The modified equation of several representative schemes is now presented. Because the 
expressions are of excessive length, only fourth-order tridiagonal schemes are considered. I ho 
equation for the ( 2-4 T ) scheme is given as 


(35) 


(36) 


(37) 


(38) 


(39) 


( 40 ) 


resulting 

modified 


dr dr d 2 r 

— — + a (\ , tt 

dt d.v d.r 2 


= “(A/)* f —T 
(j dP> 


a 

+ — 

72 


!)«'-'( A/ ) 3 


36«,.(A t) 2 - 16tf 2 (Ax) 2 (Af ) 


cl 1 /’ 

d.r { 


a 

Iso 

a 

m 

+ O f( Ax ) 1 


<M)n 2 (A /) 2 - l) 0 a 2 n r (A /)‘ 4 + 9a^(At) A + 8 ()n r /e(Ax) z (A/) 


-4<)a 2 i4 2 (Ax) 2 (A /) 2 + (Ax) 


c)x~’ 


<-) >! ' 
Ch r> 


( 41 ) 


and for t he (3-4T) HKIAV scheme as 
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t)v au d' 2 r 
W + a a7~ ap W 


a 

2 lF 


9a 2 ( A/) 3 + 16£ 2 (Ax)~(A/) 


d A U 

dx 4 


a 

+ 540 * 


+ Q 


90« 2 o v(Atf - 18a 4 (A/ ) 4 + 80o,.if 2 (Ax) 2 (A/) + 3(Ax) 4 
(Aar) 5 


cf>U 

dr> 


( 42 ] 


The term (Ax)" (At ) represents the dominant space-time error term for both the extended MaoCormack 
and RKLW schemes. For the ( 3 - 4 T) and ( 4 - 4 T) Runge-Kutt a schemes, 


dU dU d 2 V a 4 ..‘\d 4 U 

Ttr + “T7~"' T7‘- 


+ 


180 


Ma 2 a v {Aty - 6a 4 (A if + (Ax) 4 


dh’ 

dx^ 


+ 0 


(Ax)- 


(4.3) 


and 


dU dU_ d 2 U _ a 
~ + a ~ ~ ° , 'dx T ~ 360 


+ 



( 44 ) 


As expected, no space-time coupling terms exist. In each of these schemes, the first occurrence of purely 

viscous error terms (a = 0 ) is associated with -777. To retain the formal accuracy of a scheme with 

o± A) 

errors (A/)/' and (Ax) 1 /, the modified equation must only contain terms proportional to (A/)'\ (Ax)*, and 
( A/) /, (Ax) s . where r ^ pi s ^ (/, and S + R > mi n(p,q). The accuracy of the schemes is verified in 
the absence of boundaries for the convection-diffusion equation. A nonlinear viscous accuracy analysis of 
the schemes in the absence of boundaries is presented in appendix C and shows that viscous terms are 
calculated to the same accuracy as inviscid ones and that the schemes retain their advertised accuracy on 
the nonlinear problem. 


Numerical Boundary Conditions 

In each of the numerical schemes presented, a special procedure must be derived to evaluate the 
derivatives at the computational boundary points. Because accurate interior-scheme stencils are usually 
large, typically at least the derivatives at the boundary grid points require a non centered stencil. To 
preserve the formal accuracy of a spatially A'th-order accurate interior scheme on hyperbolic equations, 
the boundary and near-boundary points must be closed with stencils that are no less than ( N — 1 )th order. 
( See ref. 50 .) The procedure to derive higher order implicit and explicit boundary stencils with a symbolic 
manipulator is straightforward. (See appendix A.) Unfortunately, schemes using these higher order (fourth- 
order and greater) stencils are most often unstable (refs. Ifi and 51 ); hence, they are inappropriate to 
implement computationally. Although lower order approximations to the derivative at the boundary 
points degrade the formal accuracy of the entire numerical method, from a practical standpoint, this 
degradation Is only observed if the boundaries are a primary source of error. Therefore, if a stable, high- 
order boundary condition is available for an interior scheme, it is used. If not, the more forgiving, lower 
order formulations are used. For the viscous derivatives, viscous derivative boundary conditions in this 
study are closed to the same order as the viscous interior operator; ATh-order inviscid central derivative 
operators are closed no greater than ( A r - l)th order. An unforeseen result of this study is that closing 
the boundary points of the dissipative interior stencils cannot be done by simply using boundary closures 
derived for the centered -difference Runge-Kutt a schemes. Formally, all extended MacCormack and RKLW 
schemes are (2-2) schemes because of an interaction of the boundary and interior stencils of the inviscid 
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derivative operator. The problem is not relegated to only these two families of schemes; it affects the 
(2-4) scheme by Gottlieb and Turkel and is likely to affect other dissipative schemes. For simplicity, the 
truncation error is derived for one time step in the form of a modified equation for the seven boundary 
points used in the discretization of the equation f't + l\ v = 0 with the RKIAY integrator and the explicit 
in vise id stencil of Gottlieb and Turkel (.‘}-4R(GI )) RKIAY as follows: 


Grid point J : 


A (7 + I7A) \/-T’A( 18 + 5aA) 2 
r, + / = — Ax — (Ax) + () 


108 


162 


(Ax)" 


C-J rid point. 2: 


V, + V, = _M-5 + 7A) Aj , _ v / OA(-27+2 : 1A)( Aj ), + Q 


108 


324 


(Ax) 


Grid point A: 


( / + f' .r 


A (84 + ;>5A) \/— TA (21 + 2.1 A) 2 


1296 


-Ax + 


624 


(Ax)" + O (Ax) 


Grid point 4: 


i'i + r, = 


A (—4 + 1 7 A ) \J— 1 A (— 12 Aa A) 

Ax -| — — 


l \ .. . (-) 


Grid point 5: 


A A- 

r, + U e = Ax 

4:2 


^Ax) 2 + 0 

h I 


(Ax)‘ 


Grid point 6: 


Grid point 7: 


f I + f j — o 


(Ax 


Tin- initial condition is exp[/'(x)] with a boundary condition exp[?(— /)]. The exact solution is 
exp [?(x — /)]■ No physical boundary conditions are imposed at intermediate levels of the scheme, a 
technique which has been shown elsewhere (ref. 52) to be higher order. Because of the lack of cancellation 
at the boundary, error terms of first order are generated at the first six grid points. The RKLW scheme 
is locally first-order accurate' near the boundary and globally second-order accurate'. Use of compact 
derivative operators would spread this error over the entire domain because of the fullness of the matrix A 
instead ejf confining it to only the boundaries. 

Table 7 shows a grid refinement study of the (3-4 E) Runge-Kutta scheme given in table' 4 versus the 
(d-lk) RKIAY scheme given in table 3. Note that in this one-dimensional problem the domain contains 
only 2 full wavelengths. Degradation from the boundaries requires significant resolution: full degradation 
of the RKIAY scheme doe's not occur until resolutions on the orde r of 8000 grid points per wavelength. 
Machine precision becomes a factor as log m 1-2 < —9. 
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Time-stable wall boundary stencils for the explicit fourth-order, centered first derivative operator are 
given by 

(45) 

(46) 

(47) 

for the compact, tridiagonal fourth-order operators. Each of these equations is third-order accurate and 
results in a formally fourth-order accurate in viscid derivative operator. A stable boundary stencil for the 
third-order viscous operator is given by 

/l = ( — 1 l/l + 16/2 — 9/3 + 2/4) (48) 


/'[ = “H - 1 1/1 + 18/2 - 9 h + 2/4 ) 
h = ( -2/1 - 3/2 + 6/3 - h ) 

/I +2/2= ^(-6/1 + 4/2-f /3) 


To close the boundary point the negative complex conjugate of the Fourier image of the 

stencil at f ■ is utilized; this means for the stencil 


that 


hfi- 2 »/.//-! 27). ^/ t >2 

Aj* A.r Aj* A/ A.r 


(49) 


a— (/-l) - 


(, + l) _ Tff/ 

mj *— (/’) ^ f nr — ( /— 1 ) 

Aj* Aj* Aj* 


a I 2 ) h J n.r-(i- :i 


Aj* 


Aj* 


(50) 


where ru* is the number of grid points. Closure stencils for the viscous interior stencils must be closed 
with some care because a different number of points must be closed 011 the two sides of the comput ational 
domain. For the explicit stencils used in this study for the (N - 1) tli-order accurate viscous derivative 
operators, (A — 1) boundary points need to be closed. The forward operator requires N/’2 points to be 
closed at the right computational boundary and (N -2)/2 points at the left computational boundary. 
This rev erses when the backward viscous interior derivative operator is used. For example, the third-order 
viscous interior derivative operator is closed with expressions for / 1 , / 2 , f tlJ . on the forward stage and /j , 

fnx' fnj - 1 011 backward stage. The most interior point of these, / 2 or depending on which stage 

is being used, is closed with the negative complex conjugate of the interior stencil. Further discussion of 
boundary closures focuses on the computational domain containing the leftmost boundary points. 

St ab le , fi ft h - or de r accurate boundary stencils for the centered -difference stencils are given by 
/ 1 

h = 197/i + 690/2 - 1 380/s + 1850/4 - 1575.fr + 822/fe - 240/ 7 + 30/g ) (51) 


h - (- 18/l - 35/2 -I- 66/3 - 30/4 + 50/r, - 57/(i + 30/j - 6/s) 

h = ( + ;} /l ~ *30/*> - 20/3 + 6 O/4 - 1 0/5 + 2/o ) 

for the explicit sixth-order derivat ive operator, by 


(52) 


(53) 


/ 1 = 


1 

1680Aj* 


(-4736/1 + 14525/2 - 26250/3 + 34475 /1 - 30100/ 5 + 


16611/o- 5250/ 7 + 725/ 8 ) (54) 
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h= to! 


(- 1 :ir> / 1 - 844 h + 1635/; - 1050/1 + 575/r, - 240/, > + 69/ 7 - 10/*) 


,'20A.c 

for the compact tri diagonal sixth-order derivative, and hy 




/[ = 1 (-40:10/1 + 14790/2 - 25740/:,+ 29370/1 - 20205/r, + 7074/,; - 1250/r + 0/*) (50) 


f',= (-1416/1- 1979/2 + 2400/3 + 2820/4 - 2240/r,- 129/,; + 684/; - 200/*) (57) 

150(JA.r 

for the compact pen ta. diagonal sixth-order derivative operator. 

Results for the compact tri diagonal sixth-order derivative operator wen' first given by Carpenter 
(ref. 10). These results are referred to as 5 2 .(5‘ 2 -6-5 2 ),5' 2 schemes and are formally sixth-order accurate in 
space. By using (#/p) to designate the number of free parameters and () t 1o denote the order of accuracy 
of the first, derivative approximation to a function / at grid point 1 , we designate schemes with as 


O 


'(#//') 0 '( # //') 


— a 


niter K)r 


//.r— 1 " 


(58) 


An ( N - 1 )th-order explicit stencil representing an approximation to / requires A' grid points. By 
extending the stencil to ( A 4- 2) grid points, two degrees of freedom are added through two free parameters. 
The expression 5^ implies a fifth -order accurate stencil with two free parameters using an eight point stencil. 
To determine the minimum number of grid points needed to be closed for a centered-difference interior 
derivative operator, the band widths for matrices P and Q must be considered. If the matrix with the 
larger bandwidth has a bandwidth of then (& - l)/2 boundary grid points must be closed on each 
side of the computational domain. Where eighth-order spatial accuracy and higher is desired, additional 
bou ndar v stencils are added to allow sufficient degrees of freedom for a stable boundary closure to be 
obtained. For eighth-order accuracy, four boundary grid points were needed at each end, each with four 
degrees of freedom. Closure of tenth-order schemes may be done with six boundary grid points, each wit h 
six degrees of freedom. 

The boundaries are closed to lower order for the compact interior stencils by 


/ 1 + — 



1 

2A.r 

3 

4A7 


(-5/i +4/2 + /:*) 


(h-f l) 


(59) 


(GO) 


The explicit sixth-order interior scheme uses the same three wall points as the fourt h-order explicit scheme. 
Both schemes an' formally fourth-order accurate and are referred to as the (3,4-(5-4,3) and (3.3.4-()-4.3,3) 
schemes. 

Fifth-order accurate boundary stencils for the fifth-order viscous operator, which are used ill several 
extended Mac( ormack and RKIAV schemes, are given by 


f\ = —A (-137/1 + 300/2 - 300/;; + 200/, - 75/-, + 1 2 /,; ) (61) 

otJA.r 

/' = —A ( - 1 2/i - 65/2 + 120/;; - 00/ 4 + 20/r, - 3/,; ) ( 62 ) 

z b()A.r 

Boundary condit ions for higher order schemes are contained in appendix I). Note that stable high-order 
boundary stencils for on*' scheme are likely to be violently unstable if used with a different interior scheme. 
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All proposed central- difference stencils have been verified to have hounded left -half-plane eigenvalues for 
the matrix A. This condition is necessary for stability of Runge-Kutta schemes but does not guarantee 
stability for Navier-Stokes calculations. Table 8 gives the value of the real part of the eigenvalue of 
the matrix A. which has the largest real component for various grid sizes. In each case, the largest real 
component of any derivative operator is contained in the left half of the complex plane. This notation 
indicates that the derivative operators will be time-stable. 

In addition, the stability of all predictor-corrector (RKLW and extended MacCormack) schemes has 

also been investigated for bounded left-half-plane eigenvalues of the matrix M = ^ (A+A“ + A“A + ). 

Several of the schemes listed have eigenvalues in the right -half- plane and are unstable: for instance, all 
(2-6P), (2-8T), (3-6P), and extended MacCormack and RKLW schemes using seventh-order boundary 
stencils are unstable. The value of the largest real eigenvalue of M is given in table 9. Again, these 
schemes are only formally (2-2) schemes in the presence of numerical boundaries. 


Filters 


In reasonably well-resolved computations, numerical errors are still present and are introduced primarily 
at high wave numbers. This can readily be seen by plotting the Fourier image of the finite-difference first 
derivative versus the spectral derivative or, as is sometimes stated, modified wave number versus wave 
number. (See ref. 13.) Figure 5 shows the accuracy of various centered- difference first derivatives relative 
to the exact spectral derivative. It is immediately apparent that compact derivative operators (P ^ I) 
are more accurate than their explicit counterparts because the representation of as a polynomial is 
most accurately done as a Fade approximate. All derivative operators have no resolution at £ = tt and 
have marginal resolut ion for wave numbers near tt. Nonlinear interaction of these unresolved, lionp h vsical 
waves of various wave numbers generates higher wave-number information. When the grid is unable to 
resolve the highest wave-number information, the error is introduced into low wave numbers and eventually 
contaminates the solution. In addition, successive application of the first derivative operator to obtain a 
second-order derivat ive results in an amplification factor of unity at ^ ~ tt for centered- difference operators; 
this application facilitates what is commonly referred to as “odd-even decoupling. To suppress these 
effects, a numerical filter is used to create artificial viscosity. Several criteria exist for a useful filter. 
Eigenvalues that correspond to low wave numbers that are resolved should be virtually untouched; the 
relatively unresolved high wave numbers should be removed. Either an explicit or implicit filter can be 
chosen. Although Lele (ref. 13) uses implicit filters up to sixth order, in this study an explicit filter is 
used because 1 it is computationally more efficient and its design is more conceptually straightforward. As 
a filtering function, we seek a function that is equal to 0 at £ = tt and is equal to 1 at, £ = 0. A simple 
function to satisfy this need is 


1 — sin 


c 

2 n 

2 


(63) 


We also desire a filter whose accuracy may be chosen arbitrarily; the filtering function must have a slope 
that approaches 0 to the chosen order as £ approaches 0. Thus a (2u)th-order filter should have, to leading 
order, = £ 2n because £ tends to zero. The magnitude of the filter must also be equal to or less than 1 . 
Several discussions of filtering approaches can be found in the literature. (See refs. 19, 18, 17, and 53.) 
This discussion follows that of Eriksson (ref. 54), who derived explicit filters of second, fourth, and sixth 
order. 


For the finite-difference implementation of these filters, begin with the following definition of the explicit 
cent ral-diffore nee operator for the (2n )th-order derivative of a function /: 


/; 


(2;,) 


IZi 1- g il±l + IlzI + 

(Aa-) (2 "> (Aar)*' 2 ") 


h li+1 + ZlZ 2 , //-K3 + Z/-3 , fi+4 + fi-4 

(Aa)* 2 ") (Aa-)* 2 ") ( (Aa')< 2 ") 


(64) 


1G 


and the Fourier image given by 


* = T + 2 a cos(£) + 2 b cos(2f ) + 2c cos(3{) + 2 d cos(4£ ) + ■ * • (65) 

If consideration is given only to the second-order accurate' versions of these derivatives, the negatives of 
the coefficients are given in table 10. 

In Fourier space, the second-order accurate' stencil of | j ' ^ ~ s 2/; + is given by 

U/ = (-1)" +1 ^2 sill 0 (66) 

Table 10 shows the terms proportional to the interior stencil coefficients for filters of orders 2 (v = 1) 
t hrough 20 ( 7/ = 10) . 

If we choose' a matrix filter function D that is symmetric, then it has real eigenvalues. Because this 
filter function is based on a stencil that, when implemented with the temporal schemes discussed, has 
eigenvalues that are negative', then l r j Djj 1 1 j is always negative. This negative value guarantee's that the' 
filter is completely dissipative. The filter function is implemented on a vector U as U = ( 1 + n/jD)U. 
where U is the filtered vector; n D must be given by ( - 1 )"+ 1 2" 2 " for a (2n)th-order filter. Figure 0 shows 
the filter strength in Fourier space. 

To close the' matrix D at the boundaries and retain symmetry, skewed stencils of order v are used with 
an interior scheme of order 2n . Appendix F contains the upper left portions of the mat rix D for filters of 
orders two through twenty (v = 1 — 10). A formally tenth-order accurate scheme may now be closed with 
no concern as to how the boundary points of the filter affect the spatial operator. 

Results 

Because of the large number of schemes presented in the text, the number of permutations of boundary 
closures possible', and run times of several hours, only a few of the schemes could be run on large grids. 
The flow field of interest w as the spatially evolving, two-dimensional, compressible nitrogen-hydrogen 
shear layer. This flow field was chosen because regions of intense gradients occur both near and far 
from the flow centerline; this makes grid allocation difficult and places a heavy burden on the accuracy 
and comput ational stability of the numerical method. It is also an important prelude to t he supersonic , 
hydrogen-air. reacting shear layer found in scramjet combustors. 

Inflow conditions to all computations are Filler supersonic (ref. bo) and are described elsewhere in 
reference 1. They represent a self-similar, supersonic, nitrogen-hydrogen shear layer with \1a r — O.T), 
where Ma r is the convective Mach number. The inflow shear-layer thickness was fixed at 2 mm. A factor 
of 100 between the shear-layer thickness and the transverse' domain was used to ensure that no reflections 
from the inflow plane returned to the shearing region before the domain ended, forcing was applied at the 
inflow to the transverse velocity component in the form of a sinusoidal disturbance. Amplitudes of 4 and 
2 percent of the mean velocity were applied to the fundamental and subharmonic frequencies, respectively. 
No attempt was made to adjust other inflow variables to preserve consistency with the governing equations 
in the presence of this forcing. A result of this was a strong adjustment zone immediately down stream 
of the inflow plane where large' amplitude' disturbances propagated toward the transverse boundaries, 
potentially contaminating the shewing region with spurious boundary reflection. The forcing profile of 
the transverse velocity was in the form of a near step function centered at t he region of maximum velocity 
gradient and falling off precipitously at the shear-layer edge; this was found to minimize the amplitude of 
the disturbances impinging on the transverse boundaries. 

Principal runs were done on grids of 301 by 151 with transverse grid clustering to represent a physical 
domain of 200 by 100 mm. Only about 10 percent of the grid points was contained in t he high-shearing 
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regions of the shear layer because of this choice. The four schemes investigated were the (2-4T), (2-6T), 
(3-4T) RKLW, and (3-6T)RI\LW. No Runge-Kutta schemes are listed because several attempts were made 
to use the (4-4T) and (4-6T) schemes but they lacked sufficient dissipation to remain stable on the grid 
densities that were run. Boundary closures for the fourth- and sixth-order schemes were implemented as 
(3-4-3) and (3, 4-6-4, 3). These combinations were chosen because they were believed to represent the most 
practical schemes for these shear-layer simulations. All computations used Euler nonreflecting boundary 
conditions (refs. 36 and 57) at the two transverse boundaries and supersonic Euler boundary conditions 
at the outflow. 

Differences in accuracy between the filtered (2-4T), (2-6T), (3-4T) RKLW, and (3-6T) RKLW schemes 
were not striking. Flow-field variables like velocity, pressure, temperature, and, more importantly, the 
differentiated quantities, vorticity and dilatation, were nearly independent of the scheme. All attempted 
schemes used a tridiagonal derivative operator that was chosen to reduce the truncation error of the 
derivative operator but caused a small increase in comput ational time due to the inversion of the matrix P. 
Table 11 lists the truncation errors for various first-order derivatives, including several stencils not 
considered in any schemes in this paper. The explicit stencils have large truncation errors relative to 
the compact stencils: this is consistent with figure 5. Truncation error is minimum for the 4T, 6T, 8P, 
and 10P stencils. Further discussion of this situation is contained in appendix A. The penalty associated 
with the large CFL limits of the explicit versions of the various schemes in figures 1 through 4 is now 
clear - explicit stencils have large truncation errors. Inspection of terms that require significant resolution, 
such as V x u; and the dilatation gradient, in well-resolved simulations could help determine the efficacy 
of sixth-, eighth-, and tenth-order accurate schemes. Schemes with spatial accuracy greater than 10 may 
be readily derived by extending the current methodology; however, spectral schemes should probably be 
considered as an alternative for such highly resolved computations. 

Computational stability was found to be more sensitive to numerical method than accuracy. Extended 
Mac Cor mack schemes are more stable than the RKLW schemes and far more stable than the Runge-Kutta 
schemes. Several runs on smaller grid densities with the (4-4T) and (4-6T) schemes were completed for 
nitrogen-nitrogen shear layers; this suggests that the robustness required for disparate-mass gas mixtures 
is significantly greater. Modifying the temporally third-order Runge-Kutta. scheme with the addition of 
the predictor-corrector sequence 1o form the RKLW family of schemes noticeably increases computational 
stability. Higher order numerical boundary conditions are presented; however, they are not used because 
dissipative schemes are, by definition, lower order at the boundary. In smaller nitrogen -nitrogen simulations 
with the (4-41) and (4-6T) schemes, the formally accurate boundary closures were found to be less forgiving 
than the lower order closures and were more likely to lead to computational instability on any given grid. 
This conclusion is based on the closure response to the manner in which the shear layer was forced. 
Inadequate resolution and minimal dissipation from both the interior and boundary stencils make the 
Runge-Kutta schemes impractical for these simulations until grid densities become significantly greater 
than those chosen here. 

The appropriate choice of schemes depends strongly on the accuracy, computational stability, and OFL 
limit s of the method and the needs and resources of the user. For smaller grid densities where robustness 
was more important than CPU time, the (2-4T) scheme was very useful. The (2-6T) scheme was slightly 
slower. On larger grid sizes, such as 301 by 451, the (3-6T) RKLW scheme was preferred because it was 
slightly more accurate and actually faster than the (2-4T). A comparison of run times indicates that the 
RKLW schemes may be run to the same physical time as its corresponding extended MacCormack scheme 
(i.e.. (3-4T) RKLW versus (2-4T)), in significantly less CPU time. Relative to the (2-4T) scheme, the 
run times of the (2-6T), (3-4T) RKLW, and (3-61) RKLW schemes are 1.24,0.75, and 0.92, respectively. 
The (3-8T) RKLW scheme is likely to require only 11 percent more CPU time than the (2-4T) scheme. 
Although we did not attempt to run any of the explicit dissipative schemes, they have significantly larger 
stability envelopes, are easy to code, require no inversion of the matrix P, and are, consequently, likely to 
be very fast . Pentad iagonal schemes were presented here up to tenth order; however, these schemes may 
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not be considered competitive because of the difficulty in inverting P unt il the truncation-error penalty of 
the explicit or tridiagonal stencils is deemed sufficiently large. 

A surprising finding in this work is the effect of filters, filtering was applied to the vector U after 
all full predictor-corrector and central -difference stages in order to remove spurious information before 
it could move to lower wave numbers. The unfiltered dilatation field for the dissipative (2-11) scheme 
was badly contaminated; filtering resulted in a, noticeable improvement, figures 7(a) and f(b) show a 
section of the dilatation field from a, 101 by lb 1 simulation using the (2-4T) scheme of a nitrogen-nitrogen 
shear layer at Ma r = 0.45. with and without filtering. For comparison purposes, a. very large 401 by 601 
simulation was run with the (3-6T) RKLW scheme. (See fig. 7(c).) I he filters were found to improve the 
calculations more than any differences in temporal or spatial order of accuracy between the schemes. All 
runs were made with the tenth-order filter to avoid lowering the order of the -)*■, 5^-6-5^,5*‘ schemes. I lie 
appropriate filter order was chosen based on either the interior or boundary accuracy of the differencing 
scheme. The strongest filter that did not degrade t he accuracy of either the interior or boundary points 
was used (i.e., the interior filter order could be no less than the order of the interior scheme nor could the 
filter boundary order be less than the boundary order of the scheme). A twentieth-order filter could t hen 
be used with a. formally tenth-order scheme. 

The fact that the filters had such a. significant effect indicates that the simulations may not have been 
completely resolved. To determine whether a calculation is well resolved, a good test (in addition to 
grid refinement) is to compare filtered and unfiltered simulations. Dilatation was a particularly sensitive 
variable to resolution. The 501 by 451 calculations should be considered "model-free simulations but not 
"direct numerical simulations/ Later simulations of the nitrogen-hydrogen shear layer at M a r — 0.4o on 
the 401 by 601 grid were believed to be fully r (■'solved because contours of third derivatives of the velocity 
were not only smooth but also physically plausible. Model-free simulations that "run are no guarantee 
that all relevant scales of the problem are resolved. 

In addition, misspecification of the physical boundary conditions, wdiicli is a current topic of research, 
becomes more apparent as the accuracy of the method is increased. for sufficiently refined grids, supersonic 
Euler outflow boundary conditions an' clearly inadequate in the center of the shear layer. Dilatation 
provides a simple tool to gauge whether the non reflet* ting physical boundary conditions used on the upper 
and lower boundaries were, in fact, reflecting. Vorticity contours give virtually no indication of this 
boundary cont ami nation. 

Concluding Remarks 

An investigation was conducted of several numerical schemes that offered high spatial and temporal 
accuracv and were used in the computation of two-dimensional, spatially evolving, laminar, variable- 
density compressible shear layers. Three schemes with various temporal accuracies and arbitrary spatial 
accuracy of both the in viscid and viscous terms w r ere presented and analyzed. All integration schemes 
made use of explicit or compact finite-difference derivative operators. Extended MacCormack schemes 
retained the robustness of the original, uniformly second-order accurate method. Spatial accuracy was 
enhanced, and the stability limit was somewhat restricted. Extending the original MacCormack scheme 
resulted in longer run times; however, simulations achieved far greater spatial resolution. The (2-44 ) 
scheme, used in conjunction with a tenth-order filter, provided an accurate, computationally stable, 
general purpose numerical scheme, for large, w : e 11- re solved simulations wdiere computational stability 
(dissipation) w r as not as critical, the temporally third-order RKIAY ( R.usanov-Kut ler-Lomax-Warming) 
scheme was preferred. As with the extended MacCormack schemes, spatial accuracy of both the in viscid 
and viscous terms could be chosen freely. Stability limits of the RKLW scheme wen* large because' of 
strong resemblance to the R uuge-Kutt a cent t a.]- difference schemes. ( oinputatioual stability was achieved 
by the same space-time dissipative terms in the extended MacCormack schemes. 1 his approach made 
the RKWL schemes more stable than the Runge-Kutta schemes with the same work load per stage. An 
additional benefit of the extended MacCormack and RKLW schemes was that, computer codes written 
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with the original MaeCormack or Gottlieb-Turkel scheme may be readily upgraded to higher temporal 
and spatial accuracy with minimal effort . Third- and fourth-order Runge-Kutta schemes, although very 
accurate, possessed insufficient dissipation for the calculations conducted in this work. Filters did not add 
enough dissipation to stabilize computations with significant compressibility and variable-density effects. 
Tri diagonal difference operators were chosen for their low truncation error. Pent -adi agonal operators are 
not likely to be competitive below eighth-order accuracy. 

In each of the schemes, stability was considered for the interior operators in the convection-diffusion 
equation Uf + aU = o r U X x • Accuracy of the extended MaeCormack and RKLW schemes was verified 
for the nonlinear problem U f + F x = 0 and the viscous problem V( - [b(x)l^] x . Numerical boundary 
treatments for Runge-Kutta schemes of various orders of accuracy were chosen and evaluated to be 
asymptotically stable. Derived, formally accurate boundary conditions were given for explicit sixth-order, 
pent adi agonal sixth-order, and explicit, tridiagonal, and pentadiagonal eighth-order central -difference 
operators. Lower order closures were also presented and shown to be stable. All boundary closures 
for the extended MacC'ormack and RKLW schemes were determined to destroy the formal accuracy of the 
schemes; this problem is a serious limitation of these schemes and is likely to occur in many other common 
dissipative schemes. Apparently, this problem has gone unnoticed for over two decades. 

Damping of high wave-number, nonphysical data was accomplished for all schemes with the use of 
explicit filters, filters have been derived up to tenth order on the boundaries and twentieth order in 
the interior. These filters use explicit finite-difference stencils, are computationally efficient, and act 
predominately on high wave-number data. Results of several simulations indicate that on moderately 
well-resolved simulations, the effects of temporal and spatial accuracy differences between the schemes 
were less important than filtering effects. 


NASA Langley Research Center 
Hampton, VA 23081-2199 
July 31, 1997 
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Appendix A 

Derivation of General Stencils 

To generate the (-('titered finite-difference approximation of to Alh-order accuracy, tin' problem is 
divided into cases with k even and odd. When k = 2?) — 1 is odd. the stencil is 


xrf |2»-l) , ,(2 h- l)i , r ,(2/i-l) 

+ 4 + //- H J + 7 L //-3 


+ /,+r 1 '] + ^-2 " + //S' 


i), 


+ oifsr 11 + A>r"] + 

/i + l ~ /i-l , /;+2 ~ / i—2 //•+•) ~ //•-:) 

° (Ax)* 2 " _1 > ; (Ax)* 2 " -1 ' f (A./')* 2 "- 1 * 

, fi+4 ~ //—I . /i+5 ~ J i — - r > . i- //+(' ~ / i — (> 

(Ax)* 2 "" 1 ' ^ (Ax)< 2 " _1 > y (Ax)* 2 " -1 ' 


(Al) 


and its Fourier image is written as 

/ [* 2<7 sin(£) + 2b sin(2£) + 2r sin(H£) + 2d sin(4£) + 2r sin(r)£) -}-■••] 
[1 + 2o cos(£) + 2ii cos(2£) 4- 2-; cos(3£) + 2£ cos(4£) + • ■ •] 

Expanding the sine and cosine functions as a Taylor series gives 


* = ' £ <'■(•; 


c(2m — 1) 
( 2 jn— 1 V s 


(A2) 


( A:t) 


A -}- k' 


to 


where the functions of t ( 2 m _n are (° ^ • * * . o . /?, ■ • •). It is required t hat 4* = (*£) + 0(£ 

/// 

approximate the spectral derivative to order N. A tenth-order pentadiagona.1 approximation to/ may be 
obtained by solving six ^ ~ * ) simultaneous equations (i \ — it, = ry = t’n — iqj = 0 and i J ;i = “ 1 ) 
in six unknowns and </). Solutions do not always exist for these stencils. Similarly, when 

k — 2 it is even, then the stencil is 


■■- + *v% + o+ iuS’ + /if;;’] + + o+ a 

— T —Li— 4- •/ + 1 /i-l , J / 1’+2 + / 1 — 2 , /i+3 + //-:) 

( A.r)( 2 "> ° (A x)< 2 "> ' (Ax)* 2 "> (Ax)< 2 "> 

, /i-H + /i — 1 /;+r, + /;-r, /i-Ki + r> _ 

f (Ax)* 2 "' f (Ax)* 2 "' (Ax)* 2 ") 

and its corresponding Fourier image is given by 

[T + 2 a cos(£) + 2b cos(2£) + 2c cos(3£) + 2d cos(4£) + 2r cos(a£) T • • ■] 
* = [1 + 2 a cos(0 + 2/1 cos(20 + 2', cos(:!£) + 2* ros(-10 + • • -] 


r ( 2 ii ) 


r ( 2 II )-| 


,(2n) 


(AT, ] 


* = E (A») 

IH=l 

A tenth-order heptadiagonal approximat ion to /' 1 may he obtained hv solving eight ^ A -t ^ ) simultaneous 
equations (rp = i’-j — l'M = t’8 — t'10 = i’12 = V’U — 0 and rp = — 1) in eight unknowns (<\ . J. ~.a,b , c.d. 
and (). Again, ii is required that = (//)* + 0(^‘ V+ p. 
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For arbitrary skewed stone ils represent mg the kth derivative, we may write 


I £ r(^) IX rlk) I f(k) , r(k) , o . ,j f(k) 

i j» ( ^' ) , i* ( £ ) | <• ( k ) 

+ «!/;_ i + '»/?/;+ 1 + f; 


= T 


fi , 2 LIlz1 + “Rfi+i , b Lfi-2 + fr/?/;+-2 cl/;- 3 + <7f/;+3 

(Axl^'l (A.r)'*') (Aar)W (A./) 1 *') 


, dj/j-A + dftfi+A < ifi-h + £ff/i+5 //-/;'— 6 + Jlifi-H> 

(Aa-)t*> (Aar)*** (Ax)<*) 


A 7 


and give its Fourier image by 


= 


f [T + (o /j + «/_ ) ros({ ) + (A /j + A i ;) cos( 2£ ) + ( c/j + c j ) cos( 3£ ) + •••]') 
l + '[("/? ~ <>{.) si “(e+ ( b f{ ~ l>l.) sin(2^)+ (c ff - c f j sin(3Q+ •••] J 


f [1 + (Oft + <*/.) «w(0+ (‘ iji + t/) cos(2s ) + —]'l 
l + *[('*« ~ ( *L ) siiifO + (d/V - d/. ) sin(2£) + •-] J 


(A8) 


* = J2 V’(2;,,-2)Z (2 "' 2) + > J2 V’(2 ,„-l)£ (2 " ,_1) (A9) 

;;;=1 ;/;=l 

For the Adh derivative, vp must be either purely imaginary (/• is odd) or purely real (k is even). We now 
must solve (N + k) simultaneous equations in ( N + A’) unknowns: 17. = ( >)* (for A- even), t/’j. = (i)* -1 (for 
A- odd), and vq = 0, / = 0, i, • • • ,( A' + A 1 — 1),/ ^ A:. 

In the special case of the centered first derivative, let /> denote the number of bands in the matrix P, 
let q denote the number of bands in the matrix Q, and let ?; be A'/2. Theorder of the derivative N is then 
equal to p -f </ — 2. Simple recursion relations for the coefficients of the matrix P can be found for p < q. 
Stencils for which p > q are of marginal practical utility because they are computationally inefficient and 
have a larger truncation error than those with p = The following relations can be derived: 


(P- l)[A-(p- 1)] 

(}>+ 1) [A' - (p - 3)] 

(P-3)(P- 1 ) [A' - (/> + l)][A'-(p- 1)] 

(P + ••!)(/' + 1 ) [A r - (P ~ 5)][A' - (p - 3)] 

(p - 5)(p - 3)(p - 1) [A f - (;>+ 3)][A- - (p+ 1 )] [A' — (p— 1)] 

(P + 5 ) ( 1> + 3)(P + 1 ) [A r - (P- 7)][A - (p- 5)][A' - (p- 3)] 

(P ~ 7)(p ~ •>)(/> - 3)(p— 1) [A’ - (p + 5)][A f - {p + 3)][. ; V - {p + l)][A r - (/; - 1 )] 

(P + 7 )(p + 5)(p + 3)(// + 1 ) [A' - (p - 9)][A r - (;> - 7)] [A - - (p - 5)][A' - (p - 3)] 


(A10) 


(All) 


(A 12) 

( A 1 3 ) 


and so on. 


For example, the t hirty -second-order (A r = 32) nonadiagonal (p = 9) first derivative has the coefficients 
o = '18/65, 3 = 132/455. 7 = 176/3185. and 6 = 99/25480. 
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Similar patterns ran also be found in the matrix Q if p < q. For example. 


If we define 


then for p = 1 and 3, 


a = 


4 [yjV- (p- l) 2 ](X + 2) 
(l>+ l) 2 [A — (p — 3)p 


<I> 


= n 


r. = 




/=] 

ff 

;=i 

p- 1 
2 


(37T7 


P+ i 


b = 4>r 


[A' 1)][/V I)] 

[A' -(]>- 5)][.Y - (/>-:})] 


for p — 1,3. and o. 


. _ r . [A' - (p + 3)][A T - (P ± 1 )][A- -(/*-! )] 
f ~ 2 [A’ - ( p - 7)][jV - (p - 5)][.Y - (p - :{)] 

for p = 1 , 3. 5. and 7, 

_ [A' - (77 + 5)][.Y - ( p + 3)][A’ ~ (p + 1 )][A’ ~ (p ~ 1 )] 

' 5 [A’ - (P ~ • ) )][A’ - (p - T)][A’ - (p - 5)][.Y - {/> - :i)] 


and so on. 

If wo consider only explicit. (7; = 1. o = ;) = ■ ■■ = 0) central-difference stencils, then 

v - 0 

(l — H — 77 

n + 1 

(»-())(»- 1) 

' 2(ii+1)(h + 2) 

(»-0)(»- l)(»-2) 

3(7? + 1 )(?t + 2 )( /i + 3) 

_ (7i - ())(>? - !)(?? — 2)(n - 3) 

4 ( 71+1 ) (77 -f- 2 )( 77 + 3 )( 77 + 4 ) 

for the first derivat ive operators (first noted by Fornberg (ref. 58)). and 

7 = — 2 ( 7 / “b b c (1 -(•••• ) 

2 

(77 + 0 )( 7/ + 1 ) 


(A 14) 
(A 15) 

(Alfi) 
(A 17) 
(A 18) 

(All)) 
(A 20) 
( A 2 1 ) 
(A 22) 

(A 25) 
( A24) 
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(A25) 


2(»- 1) 

(n + 0)(n + l)(n + 2) 

2(n — 1 )(n — 2) 

(n + 0)(»> + 1 ){n + 2 )(n + 3) 

2(» — !)(»■ — 2)(» — 3) 

(n + 0 )(?? + 1 )(n + 2 ) ( 7? + 3 ){n + 4) 


( A2G) 
(A27) 


for the second derivative operators. 

Truncation error of centered first derivative operators can also he put. in general formulas. The first 
derivat ive stencil truncation error for p < q may he written as 


r>! (n - ?•)! (r)! 11/= l (2' - 1) ^(A +n 

(A’ + D! n'=i [2(»-*)+l] 


( A28) 


If we compare two different first derivative approximations that are each ATh-order accurate with the 
subscripts a and b distinguishing the two and let r a > rf r then the ratio of the truncation errors for p < q 
is given by 

(»» — »■«)! (?’(,)! n!=c, +i (2*— 1) 

("-»•/,)!(>•/,)! -o + i] (A29) 

As v (l becomes progressively larger than the ratio becomes very small; this implies, for example, that 
for a twentieth-order derivative, a. nonadiagonal derivative has far less truncation error than a tridiagonal 
derivative. Truncation error is minimized by letting p — q ( N = 4,8. 12, * • ■) or p — q — 2 (A' = 2,6. 10,* * *) 
for p < (j. This approach has been empirically verified in cases where /; < !) for all possible values of <p 
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Appendix B 

Low-Storage Runge-Kutta Methods 

Consider the equation I i + /■’. r — 0, when:' I = F(l ) and / j- = h(I j- - /Cr- 
three-stage Runge-Kutta. method may he represented in tin' equation l’i = — as 

r* = r" -« 2 ]A/''AC" 

C** = C" - «;} 1 A/" A r" - «:{-2 A/* Ah'* > 

C' ,+1 = C" - b X Xf" AC" - b 2 Xf*Ar* - b:\Xf** Al ! ** j 


The traditional 


(Bl) 


where A is tlie matrix derivative operator and A — Alternatively, in Butcher array form (ref. -12). this 


appears as 


() 


(•') 

a 2 \ 


«:il a '-V 2 


b\ b 2 Ih 


(B2) 


in some cases, storage requirements of computations should he minimized. We briefly elaborate on the 
Runge-Kutta. scheme given by W illiamson (ref. b ! ) . He derives formulas that allow storage of only two 
Valin's (2M) of the vector U, where U has a vector length M. In the low-storage format, this becomes 


= Xf"AB" - A v r* ' 
V* = 1 '" - Byq“ 
q* = Xf*AC* - A-,q" 
I'** =l’*-B. 2 q * 

q** = A/** AC** - Atf* 

C "+ 1 =U**-fhq** 


(Bit) 


or 


U* = C" - B x Xf"AV" 

V ** = V* - B 2 (\f*AV* - A 2 Xf"Ah") 
r"+‘ = C** - fl :i [\f**Af r ** - l;i( A f*AI’* - A 2 Xf"Al r ")} j 


( B4) 


Set t ing 
schemes. The 
of 


A [ to zero creates a self-starting procedure. This procedure results in a one- parameter family of 
relationship between low-storage and traditional methods may also be shown with the aid 


0 


c 2 

a 2 1 


«:*1 «:V 2 


b\ l > 2 h : > 


0 



B 1 

Ih 


Ih + ». 2 (d 2 + 1) 

[A 2 B 2 + ih } 

B> 


[A 2 (A :i lh + Ih) + Ih] 

( A 3 Ih + Ii 2 ) Ih 
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where ( ,4p ,42' ^2 • lh ) correspond to the variables (a p a 2 , 6] , b -/ . 63 ) in Williamson's work. ( See ref. 43.) 

Because the stability bounds of the two- parameter family of the third-order, three-stage Runge-Kutta 
method are independent of the parameter choice, the analysis presented in the text is valid for the low- 
storage formulas as well. 

For c*2 ^ 0, 2/3. or and ^ 0, we define 


= /(36c:J + 36c 2 - 135c 2 + 84(-2 - 12) 


-2 = 2 c 2 + 

<•2-2 



~:i = 1 2oj - 

1 8 f’2 

+ 

1 8c 2 - 

-4 = 3(>4 - 

33 <2 

+ 

1 3c 2 - 

~r, = 69c 2 - 

62 4 

+ 

28f-2 - 

-0 = 34c 2 - 

46c 2 

+ 

34c 2 - 


11^2 + 2 


The one-parameter family is given by 


B\ = c 2 

B _ 1 2C2( C -2 ~ 1 ) ( 3-2 - - 1 ) - (3-2 - - 1 ) 2 
2 ~ 144c 2 (3c 2 -2)(o 2 - 1)2 

B = -24(3c 2 -2)(c 2 -l) 2 

(3-2 - -l) 2 - T2 c 2 ( e 2 — 1 ) ( 3— 2 — - 1 ) ► 

_ ~ (bc 2 — 4c 2 + 1)-1 + 3-3 

‘ 2 “ (2c 2 + l)ci - 3 (c 2 + 2)(2c 2 - l) 2 

4 — -4- 1 + 108(2c 2 - 1)4 - 3(2r 2 - l)- 5 

‘ 3 “ 24;K2(r2- 1)4 + 72c-2-c + 72e£(2c 2 - 13) , 


(B6) 


(B7) 


wliicli can be verified, provided that none of the respective denominators vanish. Williamson s optimized 


scheme appears in Butcher form as 
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I 

3 

1 

3 

3 

3 15 

4 

~I() l(> 


1 3 8 


r; in i 




0 

1 


*1 


5 

3 

15 

( B8) 

■4 2 

B 2 

9 

To 

A:i 

B-.i 

153 

8 




~~ 128 

15 



Allen Wray of the Ames Research Center (personal communication) has considered cases where c 2 = 2/3. 
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If an extrastage' is added, a three free-parameter family of low-storage schemes may he devised. Several 
scheme's are given that are third-order accurate' for the' nonlinear problem and fourth-order accurate' for 
the linear problem {F\ r — Constant), where the stability bounds of the convection-diffusion equation are' 
the largest for four-stage schemes. If we set r-j — then 


0 

1 

i 

3 

1 

3 

5 ;i 


~V2 4 

i 

1 1 2 

i 

4 12 3 


1 0 t 


0 - — - 


3 12 1 
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3 


3 

T 

2 

1 

A 


(B9) 


If we se't — c \ . then 



Appendix C 
Nonlinear Analysis 

Consider the equation Vf + F r — [6( x)U d ] d . in a periodic domain, where b(x) is taken to be exact. 
Let B = 61 and F x = — fV The accuracy of the extended MacCorrnack (2-[2//] ) and RIvLW 

(3-[2/?]) families of schemes is now verified (the first number representing the temporal accuracy and the 
[2r? ] representing t he spatial accuracy of the overall schemes, viscous and inviscid) with viscous stencils 
whose order of accuracy is [2n-l]. This verification is done in two parts. First, the equation V f = [b{x)lF]x 
is examined to determine the accuracy of the viscous terms; second, for the equation Uf + F d . — 0 the 
nonlinear accuracy is verified. The analysis provided is somewhat different t han that given in the section 
"Formal Accuracy." Previous researchers have investigated some of these issues for spatially fourth-order 
accurate MacCorrnack schemes (ref. 9). Use of the variable n as a superscript on the vector U represents 
the nth time step; otherwise, it denotes spatial accuracy of the derivative operator. In matrix notation, 
the R.KLW scheme may be represented for the equation Uf — [b(x ) U x ].i in matrix form as 

U* = [1 + i] { \ r A~BAt}U n 
= [I + f)\ \ r A^3A~]U* 

„ (Cl) 

V* = + 

U"+ [ = Uu n + 2(1 + X r A 2 n BA 2 n )U#] 

j 


where 


A, 


7 -f^-r- If we let 

(Aj-) 


then we may define 


A + = A 2 " + X 1 
A - = A 2 " -X 
A + = A 2 " +X„ | 

A“ = A 2 " — X, J 

A 2 "BA 2 '' =21^ 
A 2 "BX ( . = Z 2 
XBA 2 '' = Z;{ 
XBX,, = Z ( 


(C2) 


(C3) 


To facilitate our analysis, note in tables 2, 4, and 7 that X r , the forward, first derivat ive matrix operator 
of (2n — l)th-order accuracy, may be rewritten as 




(04) 


where , and represent theexplicit. forward, first derivative operator 

to (2 n— l)th-order accuracy at grid point /: t he explicit , centered, first derivative operator to (2n)th-order 
accuracy at grid point i: and the explicit, centered, (2n)th derivative operator to second-order accuracy 
at grid point /, respectively. 

The matrix X may be considered by rewriting the finite- difference stencil of which it is composed 


F fi -3 + Ofi-2 + Bfi - 1 + Ofi + Bfi + 1 + Dfi+ 2 + / 7 i +3 


(C5) 
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as 


F 


+ ( ,2 %>L 2 


+ (I)+'2F) 


(°’^.), +1 + ( ,2 %,L 


+ (B + 2D + ‘AF) { [Z, ±\., V . 


(CO) 


Although we do not, utilize these results, for the special case of B = —(41) + the finite difference 
stencil in X is shown to he given by 


F 




(C7) 


where G is still given l>v G = —2(li + D + /•') or, in this particular case, (i = 2(3 /> + 8/’). By further 
enforcing D = — 0/- . equation (C7) becomes 




(C*) 


and the stencil in X then represents the second-order accurate approximation to the sixth derivative. 
The matrices A 2 ". X, and X v may he represented as 


A- = (A,) 1 2-+0 


(Ax) 


,2>i 


} 1 


X oc (Ax) 2 {-^ + O [(Ax ) 2 


X r * (Ax) 2 " 1 + ° 


(Axf 




(CD) 


The terms A^'T X. and X r that occur in both (extended Ma.c( -orniack and RKLVV schemes are now given 
bv 


Equations ((TO) imply that 


A 2 " oc (A.r) 

Xoc (Aj-) 2 
X c oc (Ax) 2 " 

Z, oc (Ax) 2 ' 

Z 2 oc (Ax ) 2,1+1 
Z, oc (Ax) 2 
Z., oc (A.r) 2,1+2 J 


(CIO) 


(Cl 1 ) 


Also note that A 2 ". X. and Xr are antisymmetric, symmetric, and symmetric, respectively. Bach 
symmetric matrix has identical diagonal elements: this implies that Z| and Zj are symmetric. Z-j and 
Z ; j are antisymmetric, and therefore 

A + BA“A“BA + = (Zi -Z.,) 2 -(Z-j - Z :! ) 2 + 2 [Z i(Z 3 - Z* ) - Z4(Z :{ + Z-j)] 1 
A+BA“ + A“BA + = 2(Z| - Z ,) J 
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The full R.KLW scheme is written as 


^ /i+l — + A |t Z| )( Z i — Z4) -f -■-X r Z\ + A f ?(I + A r Zi) 

x |(Z| -Z 4 ) 2 - (Z 2 -Z 3)' 2 + 2 [Z 1 (Z ;{ -Z 2 )-Z 4 (Z 3 + Z 2 )]|)r ,w (013) 

and the extended Mac Cor mack schemes are given as 
/ x '-2 f 

r r "+* = U+ “ {(Zl -Z 4 ) 2 -(Z 2 - Z;j ) 2 + 2 [Z|(Z;{ — Z 2 ) — Z 4 (Z;{ + Z 2 )]| 

+ a' ( .(Z| -Z 4 ))C" (014) 


With ;j[ = 1 and the RKLW schemes become 

1 / ' \ 2 1 


c 


«+i 


ib. 


:+ (a' ( .Zi) + -(a',,Zi) + -(a' ( ,z,)' - iA , l ,(i + Al.z 1 )z 4 - + A ( .Zi) 

{( Z;{ - Z 2) 2 - 2Z,Z 4 + Z'i + 2 [Z 1 ( Z :j - Z 2 ) - Z 4 (Z 3 + Z 2 )]|j U" 


and similarly, for the extended M acCormack schemes 

2 


t rn + l = (1 + (A'zi) + ^(kz { ) 2 -x c z 4 

- U, {(Z :! - Z 2 ) 2 - 2Zi Z 4 + z\ + 2 [Z [ ( Z 3 - Z-2) - Z 4 (Z 3 + Z 2 
Making use of the following relations 


)] C 


(015) 


(Cl 6) 


0" = Wtrjr,"], 

= (Aa*)“ 2 (A 2 w BA 2 ")C" 

= (Aarr 2 (Z 1 K- M (Cl 7) 


c;; = {6(r)[6(x)r.' 

= ( A X ) ■ ~ 4 ( A 2 " B A 2 " A 2 " BA 2 " )U " 

= ( Aj*) - f, (ZiZi )U" (018) 

U'ui= 

= ( Ax ) _(i ( A 2 " BA 2 " A 2 " BA 2 " A 2 " BA 2 " )U" 

= (AxJ-^ZiZjZj)!/" (Cl 9) 
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gives 

r" +1 = r" + (A/K7* + ^(A/) 2 ry; + i(A/) :, r/;, - jA '.( i + a' ( .z,)z, - 1a'; 2 (i +a',.z,) 

x |(Z ;{ - Z 2 ) 2 - 2Z1Z4 + Z'i 4 - 2 [Zi (Z 3 - Z->) - Z.| ( Z;t + Z 2 )]}r" (V 20 ) 

for the RKLW schemes and 

r" +I = e" + (A/)r 7 * + f(A/) 2 r,", - a'.z , 

- 1a' 2 {(z ;( - z 2 ) 2 - 2z 1 z , + z 2 + 2 [z,(z ;? - z.j) - z 4 (z 3 + z,)]} r" (< 21 ) 

for the extended Mac Corn lack schemes. 

The lowest order error terms present in the extended Mac( ormack scheme are those proportional 
to a' 2 (Z 3 - Z 2 ) 2 , a' 2 Zi(Z;{ - Zj), and a'.Z 4 . These throe terms have respective errors of (A/)’ 2 (Ax) 2 . 
(A/) 2 (A.r). and (A/)(A.r) 2 " (n = I [). RKLW error terms are very similar with the term (1 + A,,Zi ) being 
equal t o ( 1 + A/ ). 

Both 7 j) and Z4 arc' cross-coupling terms present m the extended Mac( ormack and RKLW schemes, 
and linear occurrences may be removed by a procedure discussed later in this appendix, luror terms that 
contain only Z\ and/or Z4 cannot be removed by this procedure. Therefore, the ( 3 - 2 E) scheme does not 
retain formal accuracy. In cases of 11 > 1 , the leading order error terms are (A/)"(A.r)" and (A/)(A.r) 

All extended MacO 'ormack and RKLW schemes except (.V 2 E) are formally accurate to their stated order 
in the absence of boundaries on the viscous problem. 

We now consider the second part of the problem, namely 

n + Fr = 0 

where F - F{V) and Fj = FfF. v = fU*. With no loss of generality, the term [fc(.r )/ V].r is absorbed into 
F x . The RKLW schemes may be represented for the equation I f — —F x as 

U* = (I — il[\f n A + )F" 

r** = (i- ;iiA/*A~)r^ 

= [r n + ihit r **-t rn )] 

u n+l = i[r"' + 2(1- A/**A 2,1 )r # ] 

whore, again. A = With 

A 2 "/A 2 " = Zr, 

A 2 "/X = Z G 
X/A 2 " = Z 7 
X/X = z„ 

Z = Zr, + Z« — Z7 — Zjs 
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the full R.KLW scheme is written as 

( ; " +1 = (i - |a/**A I 2 "/) 2 {-/?iA[(/ m + /*)A 2 " + (/"-/*)X] + tfX’rz} 

- |a/**a 2 " + ^h{-ihKif n + /*)A 2 " + (/"-f)X] 4- ,d, 2 A 2 .rz}) O" (024) 

If we expand /* and / ** in a I ay lor series with the coefficients from the Butcher array for the R.KLW 
scheme (c; = (0, then 


/* = /" + to a#)//' + — + ... 

= /" + (A i)f!' + + ... 


r = r + (r 3 Ao//* + //; + 

= /" + ^//' + ~^//; + ■ • • 


(025) 


(026) 


The expansions for /* are ident ical for the extended MacOormack and the RKLW schemes. Note that the 
error of the quantities Z 5 , Z 0 ,Z 7? and Z 8 are proportional to (Aj*) 2 , (Ax) :i . (Axf, and ( A,r) 4 . respectively. 
Inserting the values of li\ and $2 and expanding give 


'//+] 


I + ±A 2 r*A 2 M (/"+/*)A 2 '' + Ia 2 /^a 2 "(/"-/ + )X - -A ;t r*A 2 "(/"-/*)Z 

u 0 (j 


- —A/** A 2 " 
3 


-A{/" + r)A 2 " 

o 


AA (/" -DX + L\ 2 /*Z 
0 t) 


(027) 


In addition, 


o" 


= -/"o; 

= (Aj-) -1 (-/"A 2h )0" 


c//= - /,f r " + a re?), 

— ( A j- ) _ 2 ( — //' A 2 " + f"A 2n f"A 2n )U" 


(029) 


/// 


fiw'j + 2 mri^h- + rum)* - nnn^h]. 


-3, 


(-//; A 2 " + 2//'A 2,, / m A 2 " + /"A 2 " /"A 2 " - a 2 "/"A 2 "/"A 2 "H 


r" a *2/1 w/ a 2// 


a a 2// ///. A 2// 


2/1 /*/! A 2/1 A 2// \/r)| 


If we use equations (028), (029), and (030) and neglect the terms with error proportional to (A/) 4 and 
higher and cross-couple terms proportional to (A/ /(A*)*, where / + /• > 4 and higher, then after significant 
manipulation, RKLW schemes are 
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r"+ l =c" + (A/)/ 7' + + — ( a/)'V 


1 


1 


:(r Mi 


+ 

+ 

+ 


A(A/) 


A- 


r — //'X + —f"(A 2 "fX-XJA 

() 0 




A(A/h\_ A*(A/) 42m ,, |v . A 2 (A/) 


— //;x - 

r" 


-A zn f n X + 


fl l (A 2n fX-X/A 2n ) 


r" 


X'2 

—Xf'X 

(> 


(C:{i) 


The terms in brackets have an error that is proportional to (A/)’^(A,r). (A/)'*(A.r), and (A/) i (Aj*)^ ? . 
respectively: the first group must, he removed in order to retain the formal accuracy of the scheme oil the 
nonlinear problem. Ideally, all other error terms listed should also be removed, lor the linear problem, 
fj> — f n — () and A 2n fX = Xf A 2ti ; the scheme is formally (3- [2n])th-order accurate, wit ] i the highest- 

error contributed by the term ^(Xf n X )( For the nonlinear |>robleni. t he stej> from time level ( // + 1 ) 
to {v + 2) is slightly different from the step from time level (7?) to (n -f 1 ). Lot F, B, and C denob' forward, 
backward, and centered differencing, respecti vely. Tnt.il this point, consideration lias only been given to 
the HK LW scheme where t he sequence of forward-backward-centered operations a re repeated indefinitely: 
F-B-Ch • If this sequence were modified to be F-B-C. B-F-C. • * ■ . then each term for which the error 
is proportional to (A/)^(A.r) and ( A/)’ 1 (A.r) vanishes. The term ^r-(Xf n X)l n st ill remains with an error 
proportional to (A/)^( A.r)^. If the forward and backward operators are not permutated, the cross-con pling 
terms would remain and the scheme would be formally (*2-[2?/])th-order accurate. 

The full, extended Mac( 'or mack scheme is written as 

r" +[ = 1 1 - |[(r + /") a 2 " + (/"-/*)X]+ -T-z j> v" (c:vi) 


or 


" ,+I = v H + (At)r ; 1 + -(A + 


^~f,"X + y/"(A a 7X- X/A 2 ") 


r" 


+ 


\ 2 

—Xf'X 

() 


(C:W) 


Again, the two terms in the first set of brackets of equation (L33) are proportiona. 1 to (A/)"(A.r) and may 
be removed by implementing tin' scheme as an F-B. B-F, • •• sequence. Ihe error terms in tin' first set 
of brackets disappear for the linear problem and for the nonlinear problem with permutated operators, 
which leaves the last error term in the second set of brackets with an error proportional to (A/)"(A.r)^. 
In mult idimensions. this permutation would be implemented as 


F. r F y - B,B (/ B r B„ 


F,F,, F,B, y - B. t F y B,F„ - F d B„ ■ • • 


for tin' extemled MacCbrtnack schemes and 


F./F</ — B B (j 


C., C< h B./ B ij — F.rF// — C.rC iy , F,rBy — B.rF</ 


C .r C y . B.r F// — FrB /y — Cj-C */.*•• 
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for the RKLW sell cine. In the j) re sent, computations, we found it sufficient to only implement the RKLW 
scheme as 


Fj-Fy - Bj - C,C V , B x By - F,F y - C r C y , • • • 
and the extended Mac Cor mack scheme a.s 

Fj-F^ Bj-By, Bj B ( y — F ir F,y, ■ 

A discussion of multidimensional stability for finite-difference schemes can be found in the work of Beckers. 
(See ref. 59.) 

This analysis does not need to be performed for the Runge-Kutta central-difference schemes because the 
spatial and temporal accuracies are not coupled in the numerical scheme; therefore, each scheme retains 
its respective formal accuracy for the nonlinear problem. 

Further work might include replacing X with a stencil having an error that is proportional to (A.r) 4 . 
Using forward and backward differences with different weightings on all three stages of the RKLW scheme 
might minimize error terms. 
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Appendix D 

Higher Order Boundary Treatments 

When eighth-order spatially accurate, central-difference stencils are coupled with either third- or fourth- 
order temporally accurate' Runge-Kut.ta schemes, seventh-order accurate' stable' boundary stencils for the 
explicit- eight h-oreh'r accurate spatial derivative operator result anel are' given by 

/ = i05()()()At . (- 375150/1 + 1 o65172/ 2 -4(W9140/ :5 + 7180075/, 

- 8503950/ + 66767 40 / — 3378844 f- + 1024890/ 

- 157500/, + 8400/io - 840/i i + 147/u) (01) 

/ = | OMIOOA r ( 5(5400/1 ~ + 233H>45/ 3 - 4322850/, 

+ 5296900 / - 4252970 / + 2147502/; - 641320/ 

+ 94500/ - 5250/io + 045/ ] - 126/ 2 ) (02) 

/ = 2l0l)0A 7 (818(i/l “ + 181;!7 "/:t - 336385/4 

+ 407120 / - 289352/; + 114674/; - 18070/ 

- 1890/, + <530/ 1 o - 84/n + 63/ 2 ) (1)3) 


/ = (I6I8O/1 - 121338/ + 3488 10/ - 944475/ 

4 210000A.r 

+ 1234800/ - 792120/ + 323876/; - 77310/ 

+ 16800/ - 7350/, 0 + 1890/,, - 63/, 2 ) (04) 

For 1 lie t .r [diagonal eight]] -or (lor accurate derivative operator, the boundary stencils are given by 

f\ = (- 127530 /, + 738864/ - 2323230/ + 4529525/ 

1 21000A.r 

- 5668950/ + 4618740/ - 2394728/; + 730230/ 

- 105000/ + 2100/ ( | - 21l) /lJ + 189/ij, ) (D5) 

/ = •ti( ) (luAr ( 18504/1 " 060705/ - 1241625/ 

+ 1510460/- 1 18351 8/ + 578382/; - 159890/ 

+ 18900/ - 210/io+ 189/, 1 - 21/12) ( DO) 
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/ = 2IU0U0Ax (99920/l “ 818104/2 + ‘ 24{m:W h - 4643275/4 
+ 5814200/) - 4407620/; + 1947008/- - 391330/ 

- 21000/ + 18900/io - 1 050/ j + 21/ 12 ) (07) 

f a = ——(5140/ - 36240/ + 74970/ - 456225/ 

4 2 10000 Aar 2 5 4 

+ 718200/ - 474600/ + 225092/ - 65970 / 

+ 10500/ - 1050/io + 210/ 1 - 21/ia) 
and for the compact pentadiagonal eighth-order accurate derivative operator, by 

fl = 84 oooAt ' ( — 5561 52/1 + 331fi908 /2 - 10530450/ + 20529425/ 

- 25599000/ + 20734644/ - 10674356/ + 3236550/ 

- 470400/ + 12600/io - 42/ 1 + 273// 

L = (99660/ - 1065099/ + 3490410/ - 6562500/ 

2 105 000 Ax 

+ 7964320/ - 6213060/ + 301450 8/7 - 822760/ 

+ 94500/ - 945/10 + 1050/ n - 84/ 12 ) (D10) 

h = Ji 0000 a; (99920/ - 818104/ + 2402330/ - 4643275/ 

+ 5814200/ - 4407620/ + 1947008/ - 391330/ 

- 21000/ + 1 8900/ 10 - 1050/n + 21 / 2 ) (Dll) 


(D 8 ) 


(D9) 


/ = 2 ioo')OAx ( ~ 1580/1 + 17136/2 " 110670/: * “ 8662r >/4 
+ 256200/ - 101220/ + 32228/ - 6330/ 

+ 2100/ - 2100/ 10 + 1050/,i - 189/ 2 ) (D12) 


1 I 1 1 1 1 1 1 

Each of these expressions is denoted by 7 ,7 ,7 ,7 -8-7 ,7 ,7 ,7 4 schemes and is formally eighth-order 
accurate in space. We were not able to find stable, accuracy -preserving numerical boundary conditions for 
the eighth-order spatially accurate dissipative schemes. 

Seventh-order accurate boundary stencils for the seventh- order viscous operator used in t he eighth-order 
extended Mac( or mack and RKLW schemes are 

fl = 5~^(-1089/i + 2940/-4410/ + 4900/ 

- 3675/ + 1764/ -490/ + 60/) (D13) 
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/■> = -j^H-90/, - + l*50/a - 1050/, 

+ 700/5 - 315/; + «4/ 7 - 10/s) (&14) 

/•; = — (+ 10/i - 140/2 - 329/, + 700/., 

4 20 Ac 

- 350 h + 140/c - 35 /t + 4/ 8 ) (015) 

The values of f\ on the backward step and on the forward step are closed by using the seventh- 

order viscous interior stencil and the negative of its complex conjugate (in Fourier space), respectively. 
For the tenth-order spatially accurate Hunge-Kutt a schemes given in table 4, ninth-order accurate stable 
boundary stencils for the central and dissipative stencils were derived but were subject to severe ( TL 
restrict ions; hence, these are not presented. 

These boundaries may be closed to lower order also. Explicit stencils may lx* closed as (4,4. 4. c, //-!()- 
f/, j\ 4. 3,3), where x is either the fourth- or sixth-order accurate explicit stencil and // is either the fourth-, 
sixth-, or eighth-order accurate explicit stencil. Tridiagonal stencils may be implemented as (4, 4, x . //- 1 0- 
y t j-,4,3). where .r is the fourth- or sixth-order accurate tridiagonal stencil and y is either the fourth-, 
sixth-, or eighth-order accurate tridiagonal stencils. The pentad iagonal stencil is closed with (4.4..r-10- 
.r,4.4), where x is either a sixth-order accurate tridiagonal or an eighth-order accurate pen tad iagonal 
stencil. Each of these lower order closures results in a formally spatially fourth-order scheme. Other 
closures not mentioned may also be constructed. 

The viscous boundaries for the (3-10P) RKLVV scheme are closed with the following ninth-order stencils: 

f[ = ! (-7129/, + 22680/-, - 45360/-, + 70560/, - 79380/- 

1 2520 A j: 

+ 6:5504 /„ - 35280/7 + 12960/s - 2835/, + 280 /k,) (1X6) 

/.' = 1— — (- 280/, - 4329/2 + 10080/, - 11760/4 + 11760/, 

" 25 20 A. r 

— 8820 /o + 4704/- - 1680/, + 360/, - 35/ 10 ) (H17) 


/•', = - ' - -( + 35/, - 630/2 - 2754/, + 5880/, - 4410/r, 

2n2()A,r 

+ 2940/; - 1470/- + 504/„ - 105/, + 10/ 1() ) (D18) 

h = - ',- (-10/1 + 135/2 - 1080/, - 1554/, + 3780/r, 

2n20A.r 

- 1890/; + 840/t - 270/8 + 54/, - 5 /,„) (1X9) 

The value's of /- on the backward step and f\ u ._ 4 on the forward step are closed by using the ninth-order 
viscous interior stencil and the negative of its complex conjugate (in Fourier space), respectively. 
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Appendix E 

Explicit Finite- Difference Filters 

lo filter the vector U, the dissipation matrix D must be specified, including boundary points. 
Coefficients for the int erior portion have already been given in table 7 and should be familiar as elements of 
Pascals triangle. The boundary portions (upper left portion of D) of the dissipation matrix are given for 
filters of orders two through twenty in the interior and orders one through ten at the boundary. The lower 
right portion of D on a grid of A> points is given by D;j = + ^.+1— j- f° clarify the information 

presented, the full dissipation matrix is written out for the second-order filter and its first-order boundary 
points as 


+ 1 

-1 

0 

0 

0 

0 ' 


-1 

+2 

-1 

(J 

U 

0 


0 

-1 

+2 

-1 

0 

0 

(El 

0 

0 

-J 

+2 

-1 

0 

0 

0 

0 

-1 

+2 

-1 


0 

0 

0 

0 

-1 

+ 1- 



Boundary points for the second-order dissipation matrix D are 


+ 1 -1 
-1 +2 


(E2) 


The lower row and the right column are the interior operator. Similarly, the full dissipation matrix can 
be constructed for the fourth-order filter with 


for sixth order with 


for eighth order with 


for tenth order with 


-1 

+2 

-1 

+2 -5 

+4 

-1 

+4 

-6 


- + 1 

-3 

+3 

-1 ' 

-3 

+ 10 

-12 

+6 

+3 

-12 

+ 19 

-15 

.-l 

+6 

-15 

+20 - 


-1 

+4 

-6 

+4 

-1 

+4 

-17 

+28 

-22 

+8 

-6 

+ 28 

— 53 

+52 

-28 

+4 

-22 

+ 52 

-69 

+ 56 

-1 

+8 

-28 

+56 

-70 


+ 1 

-5 

+ 10 

-10 

+5 

-1 - 

—5 

+26 

— 55 

+60 

-35 

+ 10 

+ 10 

-55 

+ 126 

-155 

+ 110 

-45 

-10 

+60 

- 155 

+226 

-205 

+ 120 

+5 

-35 

+ 110 

-205 

+251 

-210 

-1 

+ 10 

-45 

+ 120 

-210 

+252. 


(E3) 


(E4) 


(E5) 


(E6) 


38 



for twelfth order with 


- -1 

+0 

- 15 

+ 20 

-15 

+6 

-1 - 

+6 

-57 

+ 96 

- 1 55 

+ 1 10 

-51 

+ 12 

-15 

+96 

-262 

+ 596 

-560 

+200 

-66 

+20 

-135 

+ 596 

-662 

+696 

-480 

+220 

-15 

+ 110 

-560 

+ 696 

-887 

+786 

-195 

+0 

-51 

+ 200 

-480 

+786 

-925 

+ 792 

. -1 

+ 12 

-66 

+ 220 

—195 

+792 

-921. 


for fourteenth order with 


+ 1 

— 7 

+ 21 

-55 

+55 

-21 

+7 

-1 ' 

—7 

+50 

-154 

+266 

-280 

+ 182 

-70 

+ 14 

+21 

-154 

+491 

-889 

+ 1001 

-721 

+529 

-91 

-55 

+266 

-889 

+ 1716 

-2 1 1 4 

+ 1756 

-966 

+ 564 

+55 

-280 

+ 1001 

-2114 

+294 1 

-2849 

+ 1981 

-1001 

-21 

+ 182 

-721 

+ 1756 

-2849 

+5582 

-2996 

+ 2002 

+7 

-70 

+ 529 

-966 

+ 1981 

-2996 

+5451 

-5005 

-1 

+ 14 

-91 

+ 564 

-1001 

+2002 

-5005 

+ 5452. 


for sixteenth order with 


' -1 

+8 

-28 

+ 56 

-70 

+56 

-28 

+8 

-1 ' 

+8 

-65 

+ 252 

-476 

+616 

—518 

+280 

-92 

+ 16 

-28 

+252 

-849 

+ 1800 

-2456 

+2184 

-1502 

+ 504 

-120 

+56 

—476 

+ 1800 

-5985 

+5720 

-5572 

+5752 

-1750 

+ 560 

-70 

+6 1 6 

-2456 

+5720 

-8885 

+9640 

-7552 

+4512 

- 1820 

+56 

—5 1 8 

+ 2184 

-5572 

+9640 

-12021 

+ 11208 

-7980 

+4568 

-28 

+280 

-1502 

+5752 

-7552 

+ 1 1208 

-12085 

+ 11452 

-8008 

+8 

-92 

+504 

- 1750 

+4512 

-7980 

+ 11452 

- 12869 

+ I 1440 

. -1 

+ 16 

-120 

+ 560 

-1820 

+4568 

-8008 

+ 11440 

- 12870. 


for eighteenth order with 


+ 1 

-9 

+56 

-84 

+ 126 

- 1 26 

+84 

-56 

+9 

-1 ' 

-9 

+82 

-555 

+ 792 

-1218 

+ 1 260 

-882 

+408 

-117 

+ 18 

+56 

-555 

+ 1578 

-5557 

+ 5528 

-5754 

+4284 

-2178 

+752 

-155 

-84 

+792 

-5557 

+ 8454 

-15941 

+ 15912 

-12810 

+7508 

-2951 

+ 816 

+ 126 

-1218 

+ 5528 

-15941 

+ 24510 

-29817 

+26496 

-17546 

+8442 

-5060 

-126 

+ 1260 

-5754 

+ 15912 

-29817 

+40186 

-4040 1 

+51052 

-18480 

+ 8568 

+84 

-882 

+4284 

-12810 

+26496 

-40401 

+47242 

-45425 

+51788 

-18564 

-56 

+408 

-2178 

+ 7508 

- 17546 

+51052 

-45425 

+48558 

-45749 

+ 51824 

+9 

-117 

+752 

-2954 

+8442 

-18480 

+5 1 788 

-45749 

+48619 

-45758 

-1 

+ 18 

-155 

+ 816 

-5060 

+8568 

-18564 

+5 1 824 

-45758 

+48620. 
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and for twentieth order with 


-1 

+ 10 

-45 

+ 120 

-210 

+ 252 

-210 

+120 

-45 

+ 10 

-1 -i 

4-10 

-101 

+400 

-1245 

+2220 

-2730 

+2352 

-1410 

+570 

-145 

+20 

— 45 

+400 

-2120 

+5860 

- 10095 

+ 13500 

-12180 

+7752 

-3435 

+ 1020 

-190 

4-1*20 

- 1245 

+5800 

-10526 

+31000 

-40935 

+38700 

-26580 

+ 13152 

-4035 

+ 1140 

-210 

+2220 

- 10095 

+31000 

-00020 

+83980 

-85035 

+03900 

-30030 

+ 15252 

-4845 

+ 252 

— 27:50 

+ 13500 

-40935 

+83980 

-124130 

+130900 

-115275 

+75300 

-38550 

+ 15504 

-210 

+2352 

-12180 

+38700 

-85035 

+ 130900 

-108230 

+ 102100 

-124725 

+77400 

-38700 

+ 120 

-1410 

+7752 

—20580 

+03900 

-115275 

+102100 

-182030 

+107500 

-125925 

+77520 

—45 

+ 570 

-3435 

+13152 

-30030 

+75300 

-124725 

+ 107500 

-184055 

+ 107950 

-125970 

+ 10 

- 145 

+ 1020 

—4035 

+ 15252 

-38550 

+77400 

-125925 

+107950 

—184755 

+ 107900 

-1 

+20 

-190 

+ 1140 

-4845 

+ 15504 

-38700 

+77520 

-125970 

+ 107900 

-184750- 


Each of these groups of boundary stencils representing ^ ^ to second-order accuracy lias a very predictable 
pattern in Fourier space with coefficients that appear in Pascal’s triangle. 
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Definitions of Symbols Used in Tables 


A.B ... . 

n.3 T.a.b. 

A a . Ab. Am 

A max 
« 

Subscripts: 

L 

It 


coefficients in predictor-corrector stencils 
coefficients in matrices P and Q 
C'FL numbers for mat rices A. B, and M 
maximum OFL number 
wave number 


left 

right 
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Table 1. Extended MacCormack Predictor-Corrector Stencils 


Scheme 

;i 

n 

mm 

i i 

c 

i) 

E 

F 

a 

Ajnax 

(2-2 E) 

0 

0 

1 

2 

1 

2 

0 

0 

0 

0 

-i 

1.00 

(2-4 E) 

0 

0 

2 

4 

0 

1 

17 

1 

7 ) 

0 

0 

2 

~7 

0.72 

(2-4 T) 

0 

1 

4 

4 

4 

£ 

A 

0 

(i 

0 

0 

2 

0 .57 

(2-6 E) 

0 

0 

4 

4 

1 

5 

4 

_ 20 

2 

7 

1 

GO 

2 

7 


0 .65 

(2-6 T) 

0 

1 

4 

7 

7 ) 

1 

4 

1 

77 

1 

45 

0 

0 


0.50 

(2-6 P) 

1 

~T74 

17 

57 

15 

U) 

5 

77 ) 

0 

0 

0 

0 

10 
_ 14 

0.48 




25 

i 

i 

7 

1 

1 

47 

0.46 

(2-8 T) 

0 

8 

42 

4 

20 

0 

•180 

2 

"77 

(2-8 P) 

1 

40 

4 

7) 

20 

27 

1 

Til 

25 

210 

1 

4 

0 

0 

7 

”7o 

0 .45 


Table 2. Skewed Viscous Stencils 


Scheme 

<h. 

H. 

h 

n L 

T 

<‘R 

b R 

( 'R 

d R 

( R 

First order forward 

0 

0 

0 

0 


1 

■ 

0 

0 

0 

Third order forward 

0 

0 

0 

1 

_ 4 


1 


0 

0 

0 

Fifth order forward 

0 

0 

1 

20 

2 

~T 

fl 

1 

2 

8 

1 

40 

0 

0 

Seventh order forward 

0 

1 

~I05 

1 

To 

4 

~5 


1 

B 

1 

1 

~777) 

0 

Ninth order forward 

1 

w7 

1 

42 

i 

7 

4 

~7> 

B 

1 


2 

71 

i 

~ 50 

1 

olid 


45 















































































Table 3. RKLW Predictor-Corrector Stencils 


Scheme 

1) 

a 

,4 

z? 

c 

D 

E 

F 

G 

^max 

(3-2E) 

0 

0 

1 

2 

1 

2 

0 

0 

0 

0 

-1 

1.59 

(3-4E) 

0 

0 

2 

3 

0 

1 

17 

2 

9 

0 

0 

4 

9 

1.34 

(3-4T) 

0 

1 

4 

3 

4 

1 

5 

0 

0 

0 

0 

2 

5 

1.06 

(3-6E) 

0 

0 

3 

4 

3 

5 

3 

~ 20 

l 

5 

1 

60 

1 

5 

2 

“5 

1.15 

(3-6T) 

0 

1 

3 

7 

0 

4 

57 

1 

36 

(j 

77 

0 

0 

164 

577 

.92 

(3-6P) 

1 

_ 114 

17 

r>7 

19 

i 

5 

0 

0 

0 

0 

2 

5 

.94 

(3-8T) 

0 


25 

32 

11 

20 

1 

20 

3 

“5 

1 

~ 480 

11 

~5o 

0 

5 

.83 

(3-8P) 

1 

30 



2 

25 

25 

216 

5 

77 

0 

0 

198 

757 

.84 

(3-1 OP) 

j 

20 

1 

2 

n 

3 

5 

101 

600 

4 

~5 

1 

600 

11 

~20 

3 

”2 

.76 
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Table 4. Fourth- and Third-Order R.unge-Kutta Stencils 


A ecu racy 








^max 

Fourth 

order 


3 

n 

a 

i, 

c 

<1 

( 

Fou It ll 
order 

■■nma 

(in;) 

(3-4E) 

0 

0 

2 

3 

1 

~ TT 

0 

0 

0 

2. Of) 

1.26 

(4-4T) 

(3-4T) 

0 

■ 

3 

4 

0 

0 

0 

0 

l .6:5 

1 .00 

('l-(iF) 

(3- OK) 

0 

1) 

3 

4 

3 

~ 20 

1 

00 

0 

0 

1 .78 

1 .09 

(4-OT) 

(3-6T) 

0 

1 

3 

7 

9 

1 

30 

0 

0 

0 

1.42 

0.87 

(4-6P) 

(TOP) 

1 

~ 11 1 

FT 

57 

F> 

19 

0 

0 

0 

0 

1 .44 

0.88 

(4-8E) 

(3-8E) 

0 

0 

1 

5 

1 

5 

4 

Too 

1 

280 

0 

1 .63 

1 .00 

d FT) 

(3-8T) 

0 

3 

8 

25 

32 

1 

20 

1 

_ 480 

0 

0 

1 .32 

0.81 

(4-8P) 

(:i- 8P) 

1 

30 

1 

9 

20 

27 

25 

210 

0 

0 

0 

1.28 

0.78 

(4-1 OK) 

(3-1 OK) 

0 

0 

5 

0 

_ 2l 

5 

84 

5 

~ 504 

i 

1200 

1 .ii3 

0.94 

(4-101) 

(3- HIT) 

0 


39 

50 

1 

Tr> 

1 

~~ 210 

1 

4200 

0 

1 .26 

0.77 

(4-10P) 

(3-10P) 

1 

20 


FT 

24 

101 

000 

1 

000 

0 

0 

1.21 

0.71 
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Table 5. Temporal Accuracy of (4-8P) Scheme on ( / + l 


( 'Oliver genet 



( 7 4 , 7 4 , 7 4 . 7 4 -8-7 4 . 7 4 . 7 4 , 7 4 ) 

CFL 

log 10 L 2 

0.8 

-5.282 

0.6 

-5.890 

0.4 

-6.038 

0.2 

-7.854 

0.1 

-9.053 

0.05 

-9.399 


Table 6. Spatial Accuracy of (4-8 P) Scheme on U( + [ 



(7 4 , 7 4 . 7 4 ,7 4 -8-7 4 ,7 4 , 7 4 . 7 4 ) 

C ’onve rgence 

Grid 

log io L -2 

rate 

41 

-4.933 


51 

-5.741 

8.34 

101 

-8.388 

8.79 

151 

-9.797 

8.00 

201 

- 1 0.000 

0.48 


Table 7. Temporal Accuracy of Third-Order Runge-Kutta and 
RKLW Schemes on U( + PT = (J 


F 

IK 

RKLW 

(3. 3-4-3, 3) 

( -divergence 

(3, 3-4-3. 3) 

( divergence 

log io L -2 

rate 

log io 1-2 

rate 

-2.512 

-2.794 

2.91 

— 1 .877 
-2.121 

2.60 

-3.081 

2.94 

-2.913 

2.02 

-4.575 

2.97 

-3.072 

2.53 

-4.804 

2.98 

-3.910 

2.43 

-5.763 

2.99 

-4.012 

2.33 

-0.005 

2.99 

-5.275 

2.20 

-7.507 

3.00 

-5.911 

2.1 1 

-8.409 

3.00 

-0.531 

2.00 

-9.309 

./ ■ 

-7.143 

2.03 

-10.099 

■■ 

-7.749 

2.01 
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Table 8. Stability of Various Boundary Closures 


Interior stencil 


max lit A ^ ) 
(7, =31) 

max h N ( A a ) 

( 7i =201) 

max ')?( A ^ ) 
( 77 = 50 1 ) 

Order 

Type 

Scheme 

Fourth 

Explicit 

(3, 3-4-3. 3) 

-2.92 x Hr"’ 

-4.32 x HP 7 

-2.75 x HP* 

Fourth 

T rid i agonal 

(3-4-3) 

-2.(50 x HP" 

-3.88 x HP 7 

-2.45 x HP* 

Sixth 

Explicit 

(3. 3. 4-6-4. 3, 3) 

-4.41 x HP" 

-0.33 x HP 7 

-3.99 x 10“* 

Sixth 

Explicit 

( 5 2 .5 2 . 5-6-5. r» 2 .. r ) 2 ) 

-1.01 x I0“ :i 

-1.17 x HP " 

-7.07 x HP 7 

Sixth 

Trid iagonal 

(3, 4-6-4. 3) 

-0.03 x nr 5 

-8.82 x IIP 7 

-5.57 x HP* 

Sixth 

Tridiagonal 

(5 2 .5 2 -0-5 2 .5 2 ) 

-1.24 x 1(P :! 

-1.50 x HP 1 ’ 

-9.00 x 10“* 

Sixth 

Pen tadi agonal 

(S.4-6-4.3) 

— 3.57 x HP" 

— 8.14 x HP 7 

-5.14 x 10“* 

Sixth 

Pen tadi agonal 


-3.30 x HP 1 

-5.08 x HP" 

-0.48 x HP (i 

Eighth 

Explicit 

(3.3.4, 4-8-4, 4, 3.3) 

- 1.22 x HP" 

- 1 .55 x HP 7 

-9.5(5 x HP 9 

Eighth 

Explicit 

(3.3,4,6-8-6.4.3.3) 

-2.18 x 10" 

-2.92 x HP 7 

-1.81 x HP* 

Eighth 

Explicit 

( 7 1 . 7 4 , 7 4 , 7 4 -8-7 4 . 7 1 . 7 4 , 7 4 ) 

-8.50 x HP 4 

-9.11 x HP (i 

-5.43 x HP 7 

Eighth 

Trid iagonal 

(3.4.4-8-4,4.31 

-5.48 x Mr" 

-7.72 x HP 7 

-4.84 x 10“ 8 

Eighth 

Trid iagonal 

(3.4 .6-8-6 , 4 .3 ) 

-7.58 x HP" 

-1.01 x 10"° 

—0.95 x HP* 

Eighth 

Trid iagonal 

(7 4 . 7 4 . 7 ’. 7 *-8-7 1 .7 1 . 7 1 . 7 1 ) 

-1.74 x 10“ :{ 

- 1.89 x HP" 

-1.13 x 10“° 

Eighth 

Pent adi agonal 

(3, 4-8-4. 3) 

-8.15 x 10- r> 

- 1.20 x l(P fi 

-7.55 x HP* 

Eighth 

Pen tadi agonal 

(7 4 ,7 4 .7 4 .7 4 -8-7 4 ,7 4 .7 4 .7 4 ) 

-4.50 x i(r ;i 

- 1.(57 x HP 4 

-8.48 x HP 0 

Tenth 

Explicit 

(3.3,4.4.4-10-4,4.4.3,3) 

-2.04 x HP (i 

-2.28 x 10“* 

-1.30 x HP 9 

Tenth 

Explicit 

(3,3.4.6.8-10-8.6,4.3.3) 

-9.58 x 1(P (I 

-1.19 x HP 7 

-7.32 x HP* 

Tent h 

Tridiagonal 

(3.4,4,4-10-4.4,4,3) 

- 1.89 x 10 -r * 

-2.38 x HP 7 

-1.4(5 x HP* 

Tenth 

Trid iagonal 

(3, 4, (5 .8- 10-8, (5.4, 3) 

— 5.78 x HP" 

-8.09 x HP 7 

-5.00 x 10“* 

Tenth 

Pen tadi agonal 

(3.4.4-10-4.4.3) 

-4.51 x HP 5 

-0.13 x HP 7 

-3.82 x HP* 

Tenth 

Pent adi agonal 

(3.4.8-10-8.4,3) 

-8.10 x 10 _r ’ 

- 1.17 x 10"° 

-7.37 x HP* 
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Table 9 . Stability of Boundary Closures for Dissipative Schemes for M — 7 (A +A7 + A-A+) 

[ 0.001 = 1 . 00 (— 3 )] 


Stencil 

(«) 

Boundary 

closure 

scheme 

Error 

max Hi A m ) 
(tt =51) 
(«) 

Error 

max ^(Ajvi) 
(// = *201 ) 
(<0 

Error 

max *tf( Am ) 
(a = 501) 
<«) 

(2-2E) [(3-2E)] 

( 1-2-1 ) 

— 3.94(— 3) [— 3.94(— 3)] 

— 2.47 ( — 4)[ — 2.47( — 4 )] 

— 3.95( — 5) [— 3.95( —5)] 

(2-4E) [(3-4E)] 

(3.3-4-3,3) 

— 3 .94 ( — 3 ) [— 3 . 94 ( —3 )] 

— 2 .4 7 ( — 4 ) [ — 2-47 ( — 4 )] 

— 3.95( —5) [— 3.90( —5)] 

( 2-4T ) [( 3-4T)] 

(3-4-3) 

— 3.94( —3) [— 3.94( —3)] 

— 2.47( — 4) [ — 2.47( — 4)] 

— 3.95(-5)[-3.95(-5)] 

( 2-6E) [(3-6E)] 

(3,3.4- G-4,3,3) 

— 3.94(— 3) [— 3.94(— 3)] 

— 2 .4 7 ( — 4)[— 2.47(— 4)] 

— 3.95( — 5) [— 3.95( —5)] 

(2-GE)[(3-GE)] 

(5 2 .5+5-G-5,5 2 ,5 2 ) 

— 3.94(— 3) [ — 3.94( — 3)] 

— 2.47( — 4) [— 2.47( — 4)] 

— 3.95( —5) [— 3.95( —5)] 

(2-GT)[(3-GT)] 

(3,4- 6-4,3) 

— 3.77(— 3) [— 3.76(— 3)] 

— 2 .44 ( — 4 ) [— 2.44(— 4)] 

— 3.93(— 5)[— 3.93(— 5)] 

(2-GT) [(3-GT)] 

(5 2 ,5 2 -(M> 2 ,5 2 ) 

-2.93( -3) [-2.93(-3)] 

— 2.27(— 4) [— 2.27(— 4)] 

— 3.81 ( —5) [— 3.81( -5)] 

(2-(iP) [(:t-GP)] 

(3.4- 6-4,3) 

+ l.«7(+2) [+2.49(4-2)] 

+ 1 ,72( +2 ) [+2. 5G( +2 )] 

+2. 12(4-2) [+3.29(4-2)] 

(2-6P) [(3-GP)] 

(5- , 5 2 ) 

+G.G5(+2)[+7.40(+2)] 

+G.G5(+2)[+7.40(+2)] 

+G G5( +2) [+7.40( +2)] 

(2-8T ) [(3-8T)] 

(3,4 t 4-M,4,3) 

+2.Gl(+0) [— 9.42(— 3)] 

+2.73(+0)[— 4,GG(— 3)] 

+2. 74( +0) [— 4.38( —3)] 

(2-8T)[(3-8T)] 

(3,4,6-8-6 ( 4,3) 

+2.Gl(+0) [— 9.42(— 3)] 

+2.73(+0) [— 4.6G(— 3)] 

+2.74( +0) [— 4.38( —3)] 

(2-8T)[(3-8T)] 

( 7h 7 b 7b 7 1 -8-7 1 , 7 1 , 7 b 7'*) 

+ 1.09( + 1)[+1.0G(+1)] 

+ 1.09( + 1) [+1.0G(+l )] 

+1.09( + 1)[+1.0G(+1)] 

(2-8P) [(3-8P)] 

(3, 4-8-4, 3) 

— 3.70(— 3) [ — 3.73( — 3)] 

— 2.43(— 4)[— 2.43(— 4)] 

— 3.92( —5) [— 3.92(— 5)] 

(2-8P) [(3-8P)] 

(7b 7b 7 1 .7 1 -8-7 l ,7 1 .7b 7 1 ) 

+ l.G8(+2)[+1.68(+2)] 

+ 1.68(+2)[+l.G8(+2)] 

+l.G8(+2)[+l.G8(+2)] 

[(2-1 OP)] 

(3, 4,4- 10-4.4,3) 

[— 3.95(— 3)] 

[ — 2.47( — 4)] 

[— 3.95( —5)] 

[(3-10P)] 

(3,4,0-10-0,4,3) 

]-3.95(— 3)] 

[ — 2.47 ( — 4 )] 

[— 3.95(— 5)] 

[(3-10P)] 

(3,4,8- 10-8,4,3) 

[— 3.95(— 3)] 

[ — 2.47 ( — 4 )] 

[-3.95(-5)] 


'‘Brackets indicate values from scheme *2: other values are from scheme 1. 
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Table JO. Interior Filter Stencil, Second- Order Accurate Stencils for 


if 2 " f 
ih :2,) 


1 



T 

a 

b 

c 

d 


/ 


■ 

■ 

J 


&} 


+2 

-i 

0 

0 

0 

0 

0 

0 

0 

0 

0 


dll' 

. d.r 1 . 


—6 

+4 

-1 

0 

0 

0 

0 

0 

0 

0 

0 


\'tf\ 


+20 

- 15 

+6 

-1 
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0 

0 

0 

0 

0 

0 


\<Ll] 


-70 

+ 56 

-28 

+ 8 

-1 

0 

0 

0 

0 

0 

0 



+25*2 

-210 

+ 120 

-45 

+10 

-1 

0 

0 

0 

0 

0 


<>' 2 l 

.777* 


-021 

+ 792 

-495 

+ 220 

-66 

+ 12 

-1 

0 

0 

0 

0 


'a 1 '/ 

777^ 


+:i4:i2 

-3003 

+2002 

-1001 

+364 

-91 

+ 14 

-1 

0 

0 

0 


'{) m f 

777™ 


-12870 

+ 1 1440 

-8008 

+4368 

-1820 

+560 

-120 

+ 16 

- 1 

0 

0 


'dllL 

jh ]b 


+48620 

-43758 

+31824 

— 1 8564 

+8568 

-3060 

+ 816 

-153 

+ 18 

-1 

0 


'o*'r 
777 > tr 


-184756 

+ 1 67960 

- 125970 

+7 7520 

-38760 

+ 15504 

— 48T) 

+ 1 140 

- 190 

+20 

-1 


Table 11. Truncation Error of First Derivative Operators 


Accuracy 

Explicit 

Tridiagonal 

1 > en t ad i agonal 

Hej)t adiagonal 

Non adiagonal 

Fourt h order 

__!_<+> 

so 4 * 

-155** 




Sixth order 

— he 

no 

'21 0(1 ^ 




Eighth order 

La* 

1 a* 

1 

20 C (. 


GOO 1 ' 

17040*' 

44 100 

220800*' 


Tenth order 

1 cl 1 

1 c 11 

1 ell 

1 ell 

•200 r)1 

“■277-/ 

124740*' 

W2120*' 

405000** 

1 4008800*' 








































Figure 1. Stability limits of extended MacCormack schemes on one- dimensional convection-diffusion 
equation as function of inviscid CLF (A) and viscous CFL (A r ) numbers. 



Viscous stability, 


Figure 2. Stability limits of HKLVV schemes on one-dimensional convection diffusion equation as function 
of inviscid CFL (A) and viscous CFL (A { ) numbers. 
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Figure 3. Stability limits of third-order Runge-K utta schemes on one-dimensional convection-diffusion 
equation as function of inviscid (TL (A) and viscous (TL (A r ) numbers. 



Figure 4. Stability limits of fourth-order Runge-K utta schemes on one-dimensional convection-diffusion 
equation as function of inviscid (TL (A) and viscous (TL (A r ) numbers. 
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Wave number, £ 

(>. Filter functions of various orders in Fourier space as funct ion of wave number 
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